1
$\begingroup$

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!

Improve this question
$\endgroup$
2
  • $\begingroup$ Please post a[i] and b[i]. $\endgroup$
    – cvgmt
    Nov 23, 2022 at 23:52
  • $\begingroup$ 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. $\endgroup$
    – Matteo
    Nov 24, 2022 at 0:28

2 Answers 2

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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 $\endgroup$
    – Matteo
    Nov 24, 2022 at 0:26
2
$\begingroup$ 
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}]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.