[HN Gopher] A Linear Algebra Trick for Computing Fibonacci Numbe...
___________________________________________________________________
A Linear Algebra Trick for Computing Fibonacci Numbers Fast
Author : todsacerdoti
Score : 51 points
Date : 2023-11-06 13:54 UTC (7 hours ago)
(HTM) web link (codeconfessions.substack.com)
(TXT) w3m dump (codeconfessions.substack.com)
| davesque wrote:
| This is also a nice way to show that the growth of fib numbers is
| exponentially bounded.
| abhi9u wrote:
| Yes, that's a good point.
| mjcohen wrote:
| Any solution to a linear recurrence relation is exponentially
| bounded by the root(s) of the characteristic polynomial.
| pfdietz wrote:
| The closed form solution as implemented will use floats with some
| fixed number of bits, right? So it cannot possible compute the
| numbers precisely except in a finite number of initial cases.
|
| 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.
| scythe wrote:
| All of this gets much easier if you notice that the smaller
| eigenvalue is less than one. So its contribution becomes
| exponentially smaller at higher F_i. Hence, you can just take
| the corresponding power of phi, multiply by some constant I
| forgot, and round to the nearest integer. No need for an exact
| formula.
| munchler wrote:
| I think you're right, and the article says "in some rare cases
| the method may produce incorrect result due to approximation
| errors".
| pfdietz wrote:
| If by "rare cases" one means "for all but a small number of
| cases".
| tedunangst wrote:
| What's "small" in this case? I only have 32GB of RAM, which
| is only enough for a vanishingly small set of Fibonacci
| numbers.
| chowells wrote:
| 79. There are only 79 Fibonacci numbers (starting at 0,
| counting the 1 twice) that are exactly represented by a
| double-precision float. That's a lot less than 32GB of
| RAM.
| __rito__ wrote:
| Yup, I implemented all these, and only naive Python array
| based approach is the most practical, Python being a GC
| language.
|
| 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).
| dietrichepp wrote:
| I think the proof for the closed-form version is accessible to
| people with a background in linear algebra. The matrix is
| diagonalizable (not _all_ matrixes are diagonalizable, but this
| one is): M = P D P^-1
|
| Here, P is some invertible matrix and D is a diagonal matrix. The
| powers of M reveal how useful this is: M^N = (P D
| P^-1) * ... * (P D P^-1)
|
| If you adjust the parentheses, you'll see N-1 terms of (P^-1 P),
| which can be removed, giving: M^N = P D^N P^-1
|
| The powers of a diagonal matrix is done by taking powers of the
| entries on the diagonal.
|
| You see the ph values (1+-[?]5)/2 in the matrix D.
|
| Diagonalization of a 2x2 matrix is simple to do on paper. The
| diagonal of D 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 D.
|
| 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.
| abhi9u wrote:
| (Author of the article here) Thanks for writing this up, it's
| much much more intuitive than what I've read everywhere else. I
| primarily looked up Knuth (TAOCP Vol 1) for this part, and he
| showed a longish proof using generating functions which looked
| too long to show in the article, and I would not have done a
| better job than him.
|
| 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.
| dietrichepp wrote:
| My own personal experience is that, in linear algebra, there
| are a lot of different ways to prove the same thing. When you
| search for a proof of something, you may find a proof that is
| overly basic and verbose, or a proof that is overly advanced
| and terse.
|
| 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))".
| abhi9u wrote:
| That's quite true. And so many books to teach it.
|
| 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. :)
| kragen wrote:
| what if you use the closed-form answer, but without any
| floating-point or other approximations
|
| 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
| abhi9u wrote:
| But that seems to be computationally same amount of work as
| the matrix form, so we get similar performance?
| kragen wrote:
| hmm, maybe? i haven't tried it so i can't be sure, but it's
| plausible
| kmill wrote:
| Yes, I believe so, up to some multiplicative factor.
|
| 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.
| tromp wrote:
| > Another interesting thing to note here is the performance of
| fib_fast and fib_closed_form algorithms. Both of these algorithms
| essentially compute the nth power of an expression and therefore
| their complexities in terms of n are same. However, the golden
| ratio based algorithm is computing the nth power of a scalar
| value which can be done much faster than computing the nth power
| of a 2X2 matrix, thus the difference in their actual run-time
| performance.
|
| 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.
| reginaldo wrote:
| For the initiated, this is SICP[1] Exercise 1.19. That exercise
| walks you through a proof without explicitly using linear
| algebra. I remember having a blast solving this back in the day.
|
| [1] https://web.mit.edu/6.001/6.037/sicp.pdf => page 61
| abhi9u wrote:
| Interesting, TIL.
|
| Thanks for sharing. :)
| abhi9u wrote:
| I just read the exercise. That's very clever.
|
| Makes me want to sit and go through the book.
| nimish wrote:
| It's easier to treat these as linear difference equations with
| constant coefficients and use the ansatz of l^{n}. This is
| morally the same thing but without the need for matrix
| manipulation.
|
| Knuth, Oren and Patashnik's Concrete Mathematics is full of these
| and is required reading for discrete math.
| SuchAnonMuchWow wrote:
| The comparison with the closed form is weird to me: since it uses
| sqrt(5), I suppose the author intended to use floating point
| numbers to compute it. But when doing so:
|
| - 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.
| abhi9u wrote:
| Author here:
|
| 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.
| SuchAnonMuchWow wrote:
| The main issue with floating point number is that it doesn't
| act like a ring
| (https://en.wikipedia.org/wiki/Ring_(mathematics)) because of
| rounding errors. For example the equality x*(y*z) = (x*y)*z
| doesn't hold because each multiplication introduce small
| rounding errors.
|
| 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.
| wrsh07 wrote:
| Please do look at the results of taking that matrix to different
| powers:
| https://www.wolframalpha.com/input?i=%5B%5B1%2C+1%5D%2C+%5B1...
|
| 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
| abhi9u wrote:
| Yes, indeed. If you see the thirty three miniatures book, it
| gives the matrix exponentiation formula for computing Fibonacci
| numbers, without any explanation behind it. To write the
| article, I tried to figure out ways it can be explained. And
| this pattern you suggested shows up nicely even when you
| express successive Fibonacci numbers in terms of the base cases
| F0 and F1. Which ultimately lead me to realize, the
| coefficients are the elements of the matrix.
| __rito__ wrote:
| Well, I am also in the reading group, and I was very happy to see
| this pop up here.
|
| 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])
|
| That's the core.
|
| 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.
| g0xA52A2A wrote:
| Reminds me of similar post with some nice visualisations
|
| https://ianthehenry.com/posts/fibonacci/
|
| https://news.ycombinator.com/item?id=36942033
| Avshalom wrote:
| A (i think) related trick:
|
| https://kukuruku.co/hub/algorithms/using-the-quick-raise-of-...
|
| https://news.ycombinator.com/item?id=8799088
| kmill wrote:
| We essentially implemented this matrix version in Lean/mathlib to
| both compute the fibonacci number and generate an efficient proof
| for the calculation.
|
| 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.
| cfen wrote:
| Why is there no mention of a generating function approach to
| derive the analytical form?
| cfen wrote:
| Why is there no mention of using a generating function approach
| to calculate the analytical function?
___________________________________________________________________
(page generated 2023-11-06 21:00 UTC)