A linear algebra trick for computing Fibonacci numbers fast
codeconfessions.substack.com
codeconfessions.substack.com
M = P Δ P^-1
M^N = (P Δ P^-1) * … * (P Δ P^-1)
M^N = P Δ^N P^-1
You see the φ values (1±√5)/2 in the matrix Δ.
Diagonalization of a 2x2 matrix is simple to do on paper. The diagonal of Δ contains the eigenvalues of M, which can be found using the quadratic formula. Proving that this is a correct way to diagonalize any diagonalizable matrix is more of a chore, but for this specific 2x2 matrix, you can just show that that you’ve found the values for P and Δ.
This is very elegant mathematically, but I would not use the closed-form solution if I wanted the exact answer, because you’d need a lot of precision, and that’s inconvenient.
The thirty three miniatures book showed a more abstract proof, which would have been accessible to people deeply familiar with the concept of basis vectors.
Finding an explanation of a linear algebra concept that strikes the right balance can be time-consuming or difficult. This is why you see, for example, a hojillion different articles explaining what quaternions are. Everybody is most comfortable at a different place on the complexity spectrum—from “a vector is (x,y,z), here are rules for how that works” to “a vector space is an abelian group V, field F, and a map in Hom(F,End(V))”.
Which is why a significant chunk of mathematics research effort and funding should go to making such explanations less time-consuming and difficult. IMHO. It may be as simple as a vetted compilation of explanations for a single idea so the individual can find the ones that work for them, or formal development of new explanations. What’s more significant: a new incremental result in some obscure corner of math, or 100M more people understanding what an eigenvalue actually means? Or a new explanation that makes a concept more relevant in an entire application area?
I am a software engineer and don't have a major in math beyond studying discrete math in college. So I really have to work hard to understand these proofs. And I certainly value proofs explained in simple terms. :)
The treatment in thirty-three miniatures describes the steps you take but doesn't mention what you are actually doing (finding eigenvalues) or leave you with any intuition for why this was a natural thing to have tried
i mean it all comes out to an integer in the end, all the radicals cancel
i see tromp already suggested a form of this
You can use a field extension here, ℚ(√5), to represent your numbers. The numbers would be represented as pairs, (a,b), which represents the number a+b√5. You can see for yourself that this set is closed under addition and multiplication.
But when you try to implement code that computes exponents with this representation, you’re probably going to fall back to the old “repeated squaring” trick. But that’s how you calculate exponents with your matrix representation anyway—and you can note that the matrix powers (for this matrix) always take the form
| a b |
| 0 a |
Sometimes these different ideas for representing the problem give us new insight or new algorithms, and sometimes it turns out to be a new perspective on the same piece of code.
You can carry out the exponentiations in the field Q[sqrt(5)], which is two-dimentional over Q. The interesting thing here is that diagonalization is trading one 2d vector space for another -- one with a very well-behaved multiplication operation (it's distributive, commutative, associative, and invertible).
There's no need to do this to quickly compute fibonacci numbers, but it's pretty neat.
Code winds up looking like:
def fib_phi(n)
((PhiRational(0,1)**n - PhiRational(1,-1)**n)/PhiRational(-1, 2)).a.to_i
end
(Technically implemented arithmetic over Q(ϕ) not Q(√5) because it simplifies the code a bit.)
- You get a constant time algorithm, since the power function is constant time on most mathematical librarie. Without entering into details on how the pow function is actually computed, on real number you can compute it with: pow(x,y) = exp(y * log(x)). It is not visible on the performance graph probably because the matrix-based algo and the closed formula are dwarfed by the linear time algo.
- You gen an incorrect result very quickly: not even considering the fact that sqrt(5) can not be represented exactly by a floating point number, the result will be out of the range of floating point numbers with n proportional to the number of bits your floating point number holds, so probably between n=50 and n=100.
At the time of writing, I was not aware that the matrix form and the closed form are same, because you can derive the golden ratio based formula from the matrix formula. Full disclaimer: I am not a math major, and I can make such mistakes :)
Although, I am curious to learn about exponentiation of real numbers. While writing this article, I was investigating how different libraries/languages implement exponentiation. While for matrix exponentiation, numpy used the algorithm I showed in the article. I also looked at the pow() function in libm of FreeBSD and it used the form you provided. I've been trying to look at an explanation behind this. Please share a pointer if you have any references.
And this kind of equations are required for the exponentiation-by-squaring method (consider for example how x**5 is computed: x**5 = (((x*x)*x)*x)*x) using a naive exponentiation approach, but also x**5 = x*(x*x)*(x*x) = x*(x*x)**2 using the exponentiation-by-squaring method).
Integer operations and matrix operations are both rings, so the exponentiation-by-squaring approach works. For floating points, it doesn't give an exact result, but could still be used in practice to get an approximate result.
The pow function doesn't use this approach however, because both arguments of the pow function are floating points numbers, and can thus compute things like pow(5, 0.5) = sqrt(5), or even 534.2**12.7. The exponentiation-by-squaring approach works only when the exponent is an integer, so we need other solutions.
How exactly the pow function (or other math function like exp, log, cos, ...) are computed is a vast research subject in itself, but if you are interested, the Handbook of Floating-Point Arithmetic (https://perso.ens-lyon.fr/jean-michel.muller/Handbook.html) is both a good introduction to the domain (and the problems associated with the finite precision of floating point number), and provides a lot of techniques to work around or minimize the rounding errors that happens when manipulating floating point numbers.
You can even make this somewhat more concrete by using the matrix representation* of a + sqrt(5) b:
[[a, 5b],
[b, a]]
*: The matrix representation is what you get if you view each value a + sqrt(5)b as a vector, and notice that multiplication by a constant is a linear map. The matrices are these linear maps, and it can be shown that matrix addition and multiplication preserve the ring structure. You can forget about the vectors at this point.
(a+b*sqrt(5)) + (c+d*sqrt(5)) = (a+c) + (b+d)*sqrt(5)
(a+b*sqrt(5)) * (c+d*sqrt(5)) = (ac+5bd) + (ad+bc)*sqrt(5)
Note that seeing the operations this way reduce a bit the number of operations compared to the matrix representation you provided: a multiplication using the formula above require 5 integer multiplications instead of 8 with the matrix [[a, 5b], [b, a]].
But this is only a constant factor improvement: the algorithm complexity is still log(n).
And the same result can be achieved by using the special structure of the matrix: for your representation, the matrix are all toeplitz matrices, and multiplying those kind of matrices require less operations than full generic matrix multiplication. For the representation in the original article, If I remember correctly, the matrices are triangular, and again you can multiply them more efficiently using this specificity.
[1] https://web.mit.edu/6.001/6.037/sicp.pdf => page 61
Makes me want to sit and go through the book.
Thanks for sharing. :)
The other approaches use arbitrary precision integers, as you've noted, so additions and multiplications are not actually fixed-cost, but scale with the number of digits you're operating on (addition scales linearly; multiplication in Python uses Karatsuba for sufficiently large N, so it scales like n^1.58).
In fact, you can't get a sublinear result with a lossless approach, since the number of digits in Fib(n) scales linearly with n, and so you need O(n) operations to just write out the digits.
This is pretty straightforward to explain with the closed-form solution: Fib(n) is roughly 1.618^n / 2.236; by taking the binary logarithm, you can see this requires about 0.7n bits to represent. (Equivalently, 0.2n decimal digits.)
Computing powers of the golden ratio should not be done as scalars, since as other posters noted that requires a lot of precision and makes for a messy computation. It's more simply done on two dimensional numbers of the form a + sqr(5)*b with a and b integers, analogous to complex numbers. Then the computational effort can be seen to equal that of the matrix powers.
https://github.com/leanprover-community/mathlib4/blob/master...
In practice this isn't very useful (the definition of Nat.fib unfolds quick enough and concrete large fibonacci numbers don't often appear in proofs) but still it shaves a bit of time off the calculation and the proof verification.
Computing the matrix power by squaring means the sizes of the integers are small until the final step, so that final step dominates the run time.
You can try to do it with arbitrary precision integers, but obviously the runtime can't be faster than the size of the answer, which technically makes it linear again.
It can be quite fast though. I especially like Julia for this. 1) because you can just tell it to use arbitrary precision ints and 2) because you can write a power function that works for everything, floats, ints, matrices with arbitrary precision integers. Which gives one of my favourite one-liners:
fib(n) = ((BigInt [1 1; 1 0])^n)[2,1]
The time will be dominated by the time to multiply integers in the final stage of the repeated squaring, which is superlinear in the number of bits. (1)
(1) for fixed word size. One might instead argue it's a function of the number of words (# of bits/word size) and that the word size must eventually increase as integers get larger, asymptotically.
It’s easy to prove that
fib(n) < 2^n
And that bound isn’t tight. A tighter one is 2log(phi) bits per step. That’s slightly less than 0.7, so that would make it less than 90 megabytes.
For all other things that I tried, I either get integer overflows (for linear algebra based implementations) or floating point errors (for the closed form).
Knuth, Oren and Patashnik's Concrete Mathematics is full of these and is required reading for discrete math.
Beautiful expo on the subject: https://eee.guc.edu.eg/Courses/Communications/COMM401%20Sign... (Chapter 10, etc.)
But, all these method, in practicality, leads to integer overflows.
Works very well for smaller Fibonacci numbers, but not for, say, the 1500th Fibonacci number.
That is very easy and practical to do instead with the naive formula using an array based approach.
fibs = [0, 1]
for i in range(2,n):
fibs.append(fibs[-2] + fibs[-1])
I have been meaning to write it up, but dealing with a bad flu.
The closed form leads to approximation errors.
I tried Julia, Numpy, Numba, Jax, Rust, but integer overflows are there in any parallelized approach.
The naive implementation with LRU caching is the fastest and most reliable so far.
I am surely missing something, you can actually use some things to overcome integer overflows, right?
I hoped to see those in this article, but they were simply not there.
Can anyone else help?
Even UInt128 is not enough for 1500th Fibonacci number. The upside of Rust, is that it doesn't fail silently. The compiler is really helpful.
Jax, Julia, Numpy confidently give wrong answers. Only way to get the feel is to eyeball a bunch of entries and check for negative numbers. Very easy to spot overflows when you are at, say, 1600th fibonacci, but one can feel confident on wrong answers on smaller numbers.
I'm not sure the maximum value matters, since I think it only actually gives you 33 decimal digits precision and the 1500th number has over 300 digits.
Although it looks like you'd need to be more careful with it too:
https://numpy.org/doc/stable/user/basics.types.html
> np.float96 and np.float128 are provided for users who want specific padding. In spite of the names, np.float96 and np.float128 provide only as much precision as np.longdouble, that is, 80 bits on most x86 machines and 64 bits in standard Windows builds.
Your code has an accidentally quadratic runtime (instead of linear). Since the array is appended to, the code regularly increases the memory region and has to move all the previous data over.
You could pre-allocate the memory as n is known ahead. Also you don't need all that memory. You only need to store the last two entries.
import sys
sys.set_int_max_str_digits(1000000)
What is this matrix doing? When you left multiply it, you're summing the elements of the left column and putting them in the top left (next fib number). Then you're moving the current number to bottom left (and top right). Incidentally, you move the previous to bottom right.
That's interesting!! It's very simple
it points out that the bit complexity of the traditional algorithm is actually quadratic. and apparently, restricted to fixed width integers, there is a constant time method.
If you work modulo M for some fixed M, then the Fibonacci sequence is periodic with period at most M^2 (because if two consecutive numbers are known, so is the whole rest of the sequence), which is a constant. Therefore, for any particular M, you can just write down the at-most-M^2 repeating sequence and what the actual period is (call that p), and then for any n the algorithm goes: reduce n mod p, and then look up the result in your big lookup table.
(Actually, this isn't constant-time, because if n is large then reducing it mod p isn't constant-time. For that matter, no algorithm can possibly be constant-time, because unless you pick a modulus for which the period is a power of 2 you always need to look at all the bits of n.)
do folks actually read knuth ? i thought his tomes were meant only for bookshelves -lol- -sarcasm-
I don't read TAOCP. I have read first few chapters of volume-1 so I am familiar with the notation and some MIX syntax.
But while writing this article, I just opened up specific topics, such as computing Fibonacci numbers in volume-1, or evaluating powers in volume-2 and I managed to understand it. Sometimes you may find back references, such as when discussing Fibonacci numbers, he references Euclid's algorithm and you may have to go back and check it out, or ignore it (depends on the context). So I believe you don't necessarily need to read from cover to the end; you can browse it based on your interest.