Using MATLAB to find a generative equation for a sequence (2024)

As @Adam Danz said, "There are several MATLAB functions here that I don't often cross paths with, so I appreciate the reminders and introductions!"

One of the things I hope to do in this sequence of discussions I plan to post, is to teach more people about some nifty things we can do In MATLAB, sometimes using tools you have never seen, or never thought of in some context, and maybe using MATLAB on problems you never thought possible.

In this post, for example, we see a tool from linear algebra like Toeplitz used to create a matrix, and then rank and null applied, all in the context of what is really a number theoretic problem. Personally, this is something I love about mathematics, that sometimes different branches of mathematics can come together in surprising ways.

Yep. In fact, it will have faster growth than any polynomial function of n, since we know the Fibonacci numbers grow exponentially, if I recall as O(phi^n/sqrt(5)) where phi comes from the golden ratio.

And any exponential function will eventually outpace any polynomial. So you know no polynomial can approximate it.

By that logic, we can even find a good approximation for the sequence, in terms of phi, valid for at least larger values of n.

I'm almost a little surprised Alpha did not try a model of that form.

Great article, @John D'Errico! There are several MATLAB functions here that I don't often cross paths with, so I appreciate the reminders and introductions!

Out of curiosity, I asked Wolfram Alpha, "What sequence is this? [1, 3, 5, 13, 25, 51, 97, 183, 337, 613]". It generated a closed-form equation, which I translated into MATLAB as follows:

I didn't bother vectorizing it, so I tested it in a loop. The results were as expected for the training set

However, when testing it beyond the initial set:

The results did not match the expected values (1101, 1959, 3457, 6059, 10557). This isn't surprising with fitted polynomials.

The obvious question is, can I prove in MATLAB, what I was able to derive from data? That is, I can assert that the recurrence relation I derived is valid for n up to 100. But is it always valid? Logically, it seems likely true, but is it ALWAYS valid?

What can I do here? First, reduce everything to Fibonacci numbers. That is, rewrite the recurrence relation into the same relation, but in terms of the Fibonacci sequence. Given the relation:

S(n) = 2*n*F(n) + 1

and the recurrence:

S(n) = S(n-5) + S(n-4) - 3*S(n-3) - S(n-2) +3*S(n-1)

First, we can simpy toss away the additive constant, since the coefficients in the recurrence relation sum to zero. If it is valid, the recurrence relation would be equally valid for

S(n) = a*n*F(n) + b

where a and b are ANY finite constants! Remember, this is a LINEAR recurrence. As such, I'll arbitrarily choose a=1, and b=0, only to make things easier to follow below.

Now rewrite the recurrence as:

-n*F(n) + 3*(n-1)*F(n-1) - (n-2)*F(n-2) - 3*(n-3)*F(n-3) + (n-4)*F(n-4) + (n-5)*F(n-5) = 0

Finally, convert this into matrix form, where I will store only the coefficients of the Fibonacci numbers. That will look like this:

Look carefully at the matrix A. You should see that rows 2 through 5 represent the basic Fibonacci recurrence. These rows merely reflect that any Fibonacci number can always be written as the sum of the previous two Fibonacci numbers.

The first row is the essence of the recurrence. Again, I've chosen to discard the additive and multiplicative constants to make it easier to read. But that changes nothing in any substantial way.

Next, comput the rank of the subarray, composed of rows 2 through 5,

It tells us no linear combination of those rows will yield a zero vector. But if we compute the rank of the entire matrix A

The rank is still 4. We learn by appending the first row, it can be represented as a sum of the other 4 rows, since the overall rank is still 4.

The conclusion is, the linear recurrence relation for S_n is valid, for any n.

I tinkered with some of the provided code. Rather than calculate a specific value in the sequence, I wrote functions to return the entire sequence up to the specified integer, n (1:n). I observed the following:

1. Using symbolic variables is about three orders of magnitude slower than using regular doubles (tested for n = 12). I am curious why you chose to use them. Given their slowness, I commented that out for subsequent testing.
2. Memoized functions appear to not work with feval, which I used in a driver function to run each variation as a test. (The feval error stated that the memoized function was unrecognized.)
3. Using a vectorized version of the function is significantly faster than the recursive Sfun. For values of n = 20, 25, and 30, the speed-up factor was about 10, 100, and 2000+, and Sfun runtimes were on the order of 0.01s, 0.1s, and 2–3s, respectively. Testing values beyond n = 30 for the recursive function was impractical. On the other hand, the vectorized function (which is included below) consistently completed on the order of 0.001s, even for n = 100. (Preallocation of Sn for larger values of n may slightly improve this function, though n is never expected to be very large.)