{"id":572,"date":"2016-05-31T20:02:29","date_gmt":"2016-05-31T20:02:29","guid":{"rendered":"https:\/\/blogs.mathworks.com\/developer\/?p=572"},"modified":"2016-06-01T22:10:09","modified_gmt":"2016-06-01T22:10:09","slug":"scaling-fibonacci","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/developer\/2016\/05\/31\/scaling-fibonacci\/","title":{"rendered":"Small Fibs Eventually Become Large Fibs – Another exercise in scaling"},"content":{"rendered":"

OK, I had way too much fun writing that last blog post<\/a> exploring algorithmic scaling. I want to explore this a bit more with a some \"real\" code, and what is more real than a function to calculate the lovely fibonacci number<\/a>, the darling of interview questions, algorithms classes, and our<\/a> own<\/a> previous<\/a> blogs<\/a>. This will be fun to explore not only to cement the method we are using to explore scalability, but also because it is fun to see how fast and elegantly the fibonacci function can be written in MATLAB.<\/p>

First, let's look at the most obvious solution based on the definition directly:<\/p>

\r\nfunction<\/span> f = fibonacci1(n)\r\n\r\nif<\/span> n < 1\r\n    f = 0;\r\nelseif<\/span> n < 3\r\n    f = 1;\r\nelse<\/span>\r\n    f = fibonacci1(n-1) + fibonacci1(n-2);\r\nend<\/span>\r\n\r\nend<\/span>\r\n\r\n<\/pre>
This is very nice and readable and reads exactly like the definition of fibonacci numbers, but we will see that unfortunately it is not very efficient at all. This is not efficient because the recursive approach needs to continually recalculate lower values of $n$      .<\/p>
Let's create a parameterized test to invoke this algorithm. As I mentioned, this solution explodes pretty quickly so I am just going to take 10 linearly spaced points from 0 to 25. Also note that the calculation is fast for small n so we need to do this in a loop to get valid measurements.<\/p>
\r\nclassdef<\/span> FirstFibonacciTest < matlab.perftest.TestCase\r\n    \r\n    \r\n    properties<\/span>(TestParameter)\r\n        size = num2cell(round(linspace(0,25,10)));\r\n    end<\/span>\r\n    \r\n    \r\n    methods<\/span>(Test)\r\n        function<\/span> testFibonacci(~, size)\r\n            for<\/span> idx = 1:1000\r\n                fibonacci1(size);\r\n            end<\/span>\r\n        end<\/span>\r\n    end<\/span>\r\nend<\/span>\r\n\r\n\r\n<\/pre>
How does it do?<\/p>
recursiveResult = runperf('FirstFibonacciTest'<\/span>);\r\n\r\nfindMins = @(resultArray) arrayfun(@(result) min(result.Samples.MeasuredTime), resultArray);\r\n\r\nax = axes('XScale'<\/span>,'log'<\/span>,'YScale'<\/span>,'log'<\/span>);\r\nylabel('Execution time (s)'<\/span>);\r\nxlabel('Problem Size'<\/span>);\r\nhold on<\/span>;\r\ngrid on<\/span>;\r\n\r\nproblemSize = linspace(0,25,10);\r\nplot(ax, problemSize, findMins(recursiveResult), 'LineWidth'<\/span>, 2)\r\nleg = legend({'Recursive'<\/span>});\r\n<\/pre>
Running FirstFibonacciTest\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n.\r\nDone FirstFibonacciTest\r\n__________\r\n\r\n<\/pre>\"\" 
$\\mathcal{O}(2^n)$. Exponential. Problem sizes larger than 25 quickly become unproductive.<\/p>
Well, what's the next approach? A common solution to calculate fibonacci numbers is to recognize that we don't need to recalculate $n-1$ and $n-2$ for every single $n$. Rather, we can just store these values into an array and calculate the number iteratively<\/i>. Note this is actually a simple form of dynamic programming as it holds on to previous solutions of the problem in memory and it uses these prior solutions to calculate subsequent values. In essence we get some pretty substantial efficency gains because we leverage memory to our advantage. How does this look? We can just store the $nth$ fibonacci number along with those for $n+1$ and $n+2$ in a three element array and work through each of the fibonacci numbers iteratively continually updating our 3 element vector as we go.<\/p>
\r\nfunction<\/span> f = fibonacci2(n)\r\n\r\nvalues = [0 1 1];\r\nfor<\/span> q = 1:n\r\n    values(1:2) = values(2:3);\r\n    values(3) = sum(values(1:2));\r\nend<\/span>\r\nf = values(1);\r\n\r\nend<\/span>\r\n\r\n<\/pre>
To observe how it scales I am going to modify our test to pass the algorithm in as another test parameter so we can reuse the test we have already written and ensure consistency of our measurements. We don't want one algorithm to be measured using a loop of 1000 and another accidently measured using a loop of 2000.<\/p>
The test looks the same with the addition of the algorithm test parameter which contains function handles to our different fibonacci implementations.<\/p>
\r\nclassdef<\/span> SecondFibonacciTest < matlab.perftest.TestCase\r\n    \r\n    \r\n    properties<\/span>(TestParameter)\r\n        algorithm = struct(...<\/span>\r\n            'Recursive'<\/span>, @fibonacci1, ...<\/span>\r\n            'Iterative'<\/span>, @fibonacci2);\r\n            \r\n        size = num2cell(round(linspace(0,25,10)));\r\n    end<\/span>\r\n    \r\n    \r\n    methods<\/span>(Test)\r\n        function<\/span> testFibonacci(~, algorithm, size)\r\n            for<\/span> idx = 1:1000\r\n                algorithm(size);\r\n            end<\/span>\r\n        end<\/span>\r\n    end<\/span>\r\nend<\/span>\r\n\r\n\r\n<\/pre>
Note that I have used the struct syntax to define the parameter which allows us to give nice names like \"Recursive\" and \"Iterative\" to each of the parameter values. This also allows us to select just these parameters and run the exact subset of interest, in this case, the iterative algorithm.<\/p>
iterativeResult = runperf('SecondFibonacciTest'<\/span>, 'ParameterName'<\/span>, 'Iterative'<\/span>);\r\nplot(ax, problemSize, findMins(iterativeResult), 'LineWidth'<\/span>, 2)\r\nleg = legend([leg.String, 'Iterative'<\/span>]);\r\n<\/pre>
Running SecondFibonacciTest\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n...\r\nDone SecondFibonacciTest\r\n__________\r\n\r\n<\/pre>\"\" 
Alright we now see the iterative approaches scales linearly ($\\mathcal{O}(n)$)) which us much nicer, but as it turns out we can do better with a closed form solution using Binet's formula<\/a>:<\/p>
$$f_n = \\frac{(1 + \\sqrt{5})^n - (1 - \\sqrt{5})^n}{2^n \\sqrt{5}}$$<\/p>
Here is how we code that up in MATLAB:<\/p>
\r\nfunction<\/span> f = fibonacci3(n)\r\nsqrt5 = sqrt(5);\r\nnum = (1 + sqrt5)^n - (1 - sqrt5)^n; \r\ndenom = (2^n)*sqrt5;\r\nf = num\/denom;\r\nend<\/span>\r\n\r\n<\/pre>
...and now we can simply add this new algorithm as new parameter to our existing test:<\/p>
properties<\/span>(TestParameter)\r\n    size = num2cell(round(linspace(0,25,10)));\r\n    algorithm = struct(...<\/span>\r\n        'Recursive'<\/span>, @fibonacci1, ...<\/span>\r\n        'Iterative'<\/span>, @fibonacci2, ...<\/span>\r\n        'ClosedForm'<\/span>, @fibonacci3);\r\nend<\/span>\r\n<\/pre>
and runperf it!<\/p>
closedFormResult = runperf('FibonacciTest'<\/span>, 'ParameterName'<\/span>, 'ClosedForm'<\/span>);\r\nplot(ax, problemSize, findMins(closedFormResult), 'LineWidth'<\/span>, 2)\r\nleg = legend([leg.String, 'Closed Form'<\/span>]);\r\n<\/pre>
Running FibonacciTest\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n.\r\nDone FibonacciTest\r\n__________\r\n\r\n<\/pre>\"\" 
That is looking very nice! In fact it demonstrates $\\mathcal{O}(1)$ constant behavior. Can't get much better than that from a complexity perspective. However, there is one more solution I'd like to show using matrix exponentials. Hmmmm, I seem to remember hearing somewhere that MATLAB is pretty good with matrices. Does that sound familiar? Well lucky us, because the solution is  quite trivial with the power of MATLAB. This is because the $nth$ fibonacci number $f_n$ is governed by the following equation:<\/p>
$$ \\left[\r\n  \\begin{array}{cc}\r\n      1 & 1 \\\\\r\n      1 & 0\r\n  \\end{array} \\right]^n\r\n  =\r\n\\left[\r\n  \\begin{array}{cc}\r\n      f_{n+1} & f_n \\\\\r\n      f_n & f_{n-1}\r\n  \\end{array} \\right] $$<\/p>
Yes, that's right, the solution in MATLAB is just that simple.<\/p>
\r\nfunction<\/span> f = fibonacci4(n)\r\n\r\nF = [1 1;1 0]^n;\r\nf = F(1,2);\r\n\r\nend<\/span>\r\n\r\n<\/pre>
How does this scale? Add the test and run it.<\/p>
properties<\/span>(TestParameter)\r\n    size = num2cell(round(linspace(0,25,10)));\r\n    algorithm = struct(...<\/span>\r\n        'Recursive'<\/span>, @fibonacci1, ...<\/span>\r\n        'Iterative'<\/span>, @fibonacci2, ...<\/span>\r\n        'ClosedForm'<\/span>, @fibonacci3, ...<\/span>\r\n        'MatrixExponential'<\/span>, @fibonacci4);\r\nend<\/span>\r\n<\/pre>
matrixResult = runperf('FibonacciTest'<\/span>, 'ParameterName'<\/span>, 'MatrixExponential'<\/span>);\r\nplot(ax, problemSize, findMins(matrixResult), 'LineWidth'<\/span>, 2)\r\nleg = legend([leg.String, 'Matrix Exponential'<\/span>]);\r\n<\/pre>
Running FibonacciTest\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\n..........\r\nDone FibonacciTest\r\n__________\r\n\r\n<\/pre>\"\" 
Look at that! The closed form and the matrix solutions both show very scalable behavior, but the matrix form is not quite as good, showing logarithmic ($\\mathcal{O}(\\log n)$) scaling instead of constant and the constant factor is greater for the matrix case. However, there is another benefit enjoyed by the matrix exponential solution but not by the closed form solution. To show this, I am going to add a quick spot check to this test to confirm that we are actually getting the right answer. Note for this post I am just going to add another test method which validates a single size, but this could be done in the existing test method across all sizes if desired.<\/p>
\r\nclassdef<\/span> VerifiedFibonacciTest < matlab.perftest.TestCase\r\n    \r\n    \r\n    properties<\/span>(TestParameter)\r\n        size = num2cell(round(linspace(0,25,10)));\r\n        algorithm = struct(...<\/span>\r\n            'Recursive'<\/span>, @fibonacci1, ...<\/span>\r\n            'Iterative'<\/span>, @fibonacci2, ...<\/span>\r\n            'ClosedForm'<\/span>, @fibonacci3, ...<\/span>\r\n            'MatrixExponential'<\/span>, @fibonacci4);           \r\n    end<\/span>\r\n    \r\n    \r\n    methods<\/span>(Test)\r\n        function<\/span> testFibonacci(~, algorithm, size)\r\n            for<\/span> idx = 1:1000\r\n                algorithm(size);\r\n            end<\/span>\r\n        end<\/span>\r\n        function<\/span> testCorrectAnswer(testCase, algorithm)\r\n           testCase.verifyEqual(algorithm(16), 987, ...<\/span>\r\n               'The 16th fibonacci number should be 987'<\/span>); \r\n        end<\/span>\r\n    end<\/span>\r\nend<\/span>\r\n\r\n\r\n<\/pre>
Let's run it with runtests<\/tt><\/b> to confirm correctness.<\/p>
result = runtests('VerifiedFibonacciTest'<\/span>)\r\n<\/pre>
Running VerifiedFibonacciTest\r\n..........\r\n..........\r\n..........\r\n..........\r\n..\r\n================================================================================\r\nVerification failed in VerifiedFibonacciTest\/testCorrectAnswer(algorithm=ClosedForm).\r\n\r\n    ----------------\r\n    Test Diagnostic:\r\n    ----------------\r\n    The 16th fibonacci number should be 987\r\n\r\n    ---------------------\r\n    Framework Diagnostic:\r\n    ---------------------\r\n    verifyEqual failed.\r\n    --> The values are not equal using \"isequaln\".\r\n    --> Failure table:\r\n                Actual    Expected           Error               RelativeError    \r\n                ______    ________    ____________________    ____________________\r\n            \r\n                987       987         4.54747350886464e-13    4.60736930989325e-16\r\n    \r\n    Actual double:\r\n             9.870000000000005e+02\r\n    Expected double:\r\n           987\r\n\r\n    ------------------\r\n    Stack Information:\r\n    ------------------\r\n    In \/mathworks\/inside\/files\/dev\/tools\/tia\/mdd\/performanceScaling\/VerifiedFibonacciTest.m (VerifiedFibonacciTest.testCorrectAnswer) at 21\r\n================================================================================\r\n..\r\nDone VerifiedFibonacciTest\r\n__________\r\n\r\nFailure Summary:\r\n\r\n     Name                                                           Failed  Incomplete  Reason(s)\r\n    ============================================================================================================\r\n     VerifiedFibonacciTest\/testCorrectAnswer(algorithm=ClosedForm)    X                 Failed by verification.\r\n    \r\n\r\nresult = \r\n\r\n  1x44 TestResult array with properties:\r\n\r\n    Name\r\n    Passed\r\n    Failed\r\n    Incomplete\r\n    Duration\r\n    Details\r\n\r\nTotals:\r\n   43 Passed, 1 Failed, 0 Incomplete.\r\n   7.3549 seconds testing time.\r\n\r\n<\/pre>
...and there is the problem with the closed form solution. It is subject to floating point precision errors. Of course we can simply make this test pass by leveraging a tolerance in verifyEqual<\/p>
testCase.verifyEqual(algorithm(16), 987, 'RelTol'<\/span>, 1e-14...<\/span>\r\n    'The 16th fibonacci number should be 987'<\/span>);\r\n<\/pre>
but the fact remains that using the closed form solution we compute something that is close to the $nth$ fibonacci number, but the matrix solution computes it exactly. Which one to use depends on your application, but given the constant factor is so small and logarithmic behavior approaches constant for large problem sizes I think I would tend to choose the exact solution. How about you?<\/p>