What are some good ways to compute terms in recurrence relations?

General Usage
1

Doing a google search on “recurrence relations Julia” shows up MathToolkit.jl which has not been updated for a while and it’s got a neat function for detecting recurrence relations but not for solving them given the coefficients.

For recurrence relations that “reduces”, Memoize.jl works really well. E.g. Fibonaci sequence is just

using Memoize
@memoize function fib(n)
    if n <= 2
        return 1
    else
       return fib(n-1)+fib(n-2)
    end
end

And that’s it! However, I have issues with recurrence relations that don’t “reduce” e.g. if f(n) depends on f(m) where m > n.

E.g. say I am looking for the expected number of steps needed to reach finished line and the steps taken are controlled by a die where rolling 1 means take one step back and rolling 6 means go forward one step and rolling anything else means staying put. Then I get a recurrence relations like this

f(n) = 1 + f(n-1)//6 + f(n+1)//6 + f(n)*4//6

where n is the number of steps still to go in order to reach finish line.

I am just wondering if there’s a good way to solve this in Julia simply. I wrote some recursive function to solve it.

2

This isn’t a general answer, but it seems like in the example there you could just solve for f(n+1) and then have a recurrence relation only in terms of lower values, at which point you can just memoize.

3

Why didn’t I think of that? The actual problem I am solving is this #227 The Chase - Project Euler

So wondering if that applies there.

4

I don’t think you can do that in this case, since you don’t have enough boundary conditions.

I would be interested to see that function.

What I understand is that you allow positions n=1, \ldots, L, where L is the finish line, and I presume that if you are at n=1 and you try to go left, you stay at n=1 (i.e. you “bounce back” from a wall); in this way there are only a finite number of allowed positions n.

Your f(n) is the expected or mean time to hit the finish line L, starting from site n; this is a hitting time or (mean) first-passage time. In particular, f(L) = 0, since you’re already there.

If you write down all the equations for all values of n, you have a system of linear equations for the unknowns f(1), f(2), …, which you can solve using the \ operator.

A more general technique, in particular, for problems in which there are an infinite number of possible values of n, is to use a generating function, i.e. define

F(z) := \sum_n f(n) z^n

The recurrence relation, together with the boundary conditions, gives you an algebraic equation for F(z), which you can solve. You can then (sometimes) extract the results f(n) as the coefficient of z^n in the resulting expression.

2 Likes
5

Check it out. It solves the more complicated version not the simplified one I have in my original post.

# this is an alternate solution that doesn't rely on matrices so is more memory efficient

init_rhs() = begin
    res = Rational{BigInt}[0//1 for i in 1:50]
    res
end

function update227!(rhs_coef, const_part, n, multiplier)
    new_const_part, add_rhs_coef = multiplier .* chase(n)
    rhs_coef .+= add_rhs_coef

    const_part + new_const_part
end

using Memoize

@memoize function chase(dist)
    # C*d_{dist} = A*d_{dist} + B*d_{dist-1} + D*d_{dist-2}
    # lhs_coef = C - A
    lhs_coef = 1//1
    const_part = 1//1

    rhs_coef = init_rhs()

    if dist == 50
        rhs_coef[48] += 2//36
        rhs_coef[49] += 16//36
        rhs_coef[50] += 1//2
    elseif dist == 49
        rhs_coef[47] += 1//36
        rhs_coef[48] += 8//36
        rhs_coef[49] += 1//2 + 1//36
        rhs_coef[50] += 8//36
    elseif dist == 1
        rhs_coef[1] += 1//2 + 1//36
        rhs_coef[2] += 8//36
        rhs_coef[3] += 1//36
    elseif dist == 2
        rhs_coef[1] += 8//36
        rhs_coef[2] += 1//2
        rhs_coef[3] += 8//36
        rhs_coef[4] += 1//36
    else
        rhs_coef[dist-2] += 1//36
        rhs_coef[dist-1] += 8//36
        rhs_coef[dist] += 1//2
        rhs_coef[dist+1] += 8//36
        rhs_coef[dist+2] += 1//36
    end

    if dist != 1
        while any(!=(0), @view rhs_coef[1:dist-1])
            for i in 1:dist-1
                if rhs_coef[i] != 0
                    mult = rhs_coef[i]
                    rhs_coef[i] = 0
                    const_part = update227!(rhs_coef, const_part, i, mult)
                end
            end
        end
    end

    lhs_coef -= rhs_coef[dist]
    rhs_coef[dist] = 0//1

    const_part / lhs_coef, rhs_coef ./ lhs_coef
end
1 Like
6

I am well aware of generating functions from my discrete math course. Just wondering if there’s a way to input the recurrence relations coefs and a program can solve it for me using whichever method including generating functions.

7

This was my original idea. However, it fails on large matrices on resource constraint machine (like those on HackerRank) and I just realised that you don’t need a matrix for this more specialised problem, hence the question. You should be able to solve this in a more RAM-efficient way. As I have shown above.

8

For this particular problem the matrix is symmetric and tridiagonal. Using the SymTriDiagonal type should be very efficient (in memory and time) even for large systems.

1 Like
9

That’s an interesting question (which I don’t have a good answer to). Maybe you could use ModelingToolkit somehow to get the structure of the recurrence relation.

The book Generatingfunctionology may be useful.

10

These are very elegant techniques but can run into practical problems when implemented numerically (eg a very nice or even perfectly accurate calculation turning into something less well-conditioned).

Of course these problems can be mitigated, but for small n I would just use the recurrence relation, maybe implemented with a loop, or memoized. Or for non-triangular cases, use a linear system as you suggested.

1 Like