One of the main features of Octave we will discuss is vectorization. To understand it, we need some background material on Linear Algebra. If you've taken a linear algebra course recently, this should be easy, otherwise you will probably need to review
Vector Operations, particularly
Matrix Multiplication Review
We will mainly rely on the Octave Interpreter Reference. A more tutorial style guide is the Gnu Octave Beginner's Guide, which is available in ebook form from the UNB library.
We'll refer to the online text by Robert Beezer for linear algebra background. We can happily ignore the stuff about complex numbers.
There is a GUI accessible from Applications -> FCS -> GNU Octave, or by running from the command line
% octave --gui
There is also a REPL accessible from the command line by running
% octave
To sanity check your octave setup, run the following plots
>> surf(peaks)
>> contourf(peaks)
Here is a JavaScript recursive function for Fibonacci:
function fib(n) {
if (n<=0)
return 0;
if (n<=2)
return 1;
else
return fib(n-1)+fib(n-2);
}
Let's translate this line by line into an Octave function.
Save the following in ~/cs2613/labs/L20/recfib.m; the name of the
file must match the name of the function.
function ret = recfib(n)
if (n <= 0)
ret = 0;
elseif (n <= 2)
ret = 1;
else
ret = recfib(n-1) + recfib(n-2);
endif
endfunction
Like the other programming languages we covered this term, there is a built in unit-test facility that we will use. Add the following to your function
%!assert (recfib(0) == 0);
%!assert (recfib(1) == 1);
%!assert (recfib(2) == 1);
%!assert (recfib(3) == 2);
%!assert (recfib(4) == 3);
%!assert (recfib(5) == 5);
Note the %! assert. These are unit tests that can be run with
>> test fib
The syntax for %!assert is a bit fussy, in particular the
parentheses are needed around the logical test.
A powerful technique for making recursive A more problem specific approach (sometimes called dynamic programming) is to fill in values in a table.
Save the following in ~/cs2613/labs/L20/tabfib.m. Complete the missing line by
comparing with the recursive version, and thinking about the array indexing.
function ret = tabfib(n)
table = [0,1];
for i = 3:(n+1)
table(i)=
endfor
ret = table(n+1);
endfunction
%!assert (tabfib(0) == 0);
%!assert (tabfib(1) == 1);
%!assert (tabfib(2) == 1);
%!assert (tabfib(3) == 2);
%!assert (tabfib(4) == 3);
%!assert (tabfib(5) == 5);
Let's measure how much of a speedup we get by using a table.
Of course, the first rule of performance tuning is to carefully test
any proposed improvement. The following code gives an extensible way
to run simple timing tests, in a manner analogous to the Python
timeit method, whose name it borrows.
# Based on an example from the Julia microbenchmark suite.
function timeit(func, argument, reps)
times = zeros(reps, 1);
for i=1:reps
tic(); func(argument); times(i) = toc();
end
times = sort(times);
fprintf ('%s\tmedian=%.3fms mean=%.3fms total=%.3fms\n',func2str(func), median(times)*1000,
mean(times)*1000, sum(times)*1000);
endfunction
We can either use timeit from the octave command line, or build a little utility function like
function bench
timeit(@recfib, 25, 10)
timeit(@tabfib, 25, 10)
endfunction
What are the new features of Octave used in this sample code?
Roughly how much speedup is there for using tables to compute the 25th fibonacci number?
As a general principle we may want to reduce the amount of storage
used so that it doesn't depend on the input. This way we are not
vulnerable to e.g. a bad input crashing our interpreter by running out
of memory. We can improve tabfib by observing that we never look
more than 2 back in the table. Fill in each blank in the following
with one variable (remember, ret is a variable too).
function ret = abfib(n)
a = 0;
b = 1;
for i=0:n
ret = _ ;
a = _ ;
b = _ + b ;
endfor
endfunction
%!assert (abfib(0) == 0);
%!assert (abfib(1) == 1);
%!assert (abfib(2) == 1);
%!assert (abfib(3) == 2);
%!assert (abfib(4) == 3);
%!assert (abfib(5) == 5);
n==0 (first blank)n==1 (second blank)i, b is assigned to (i+2)nd Fibonacci number.function bench2
timeit(@tabfib, 42, 10000)
timeit(@abfib, 42, 10000)
endfunction
The following is a
well known
identity about the Fibonacci numbers F(i).
[ 1, 1;
1, 0 ]^n = [ F(n+1), F(n);
F(n), F(n-1) ]
This can be proven by induction; the key step is to consider the matrix product
[ F(n+1), F(n); F(n), F(n-1) ] * [ 1, 1; 1, 0 ]
Since matrix exponentiation is built-in to octave, this is particularly easy to implement the formula above in octave
Save the following as ~/cs2613/labs/L20/matfib.m, fill in the
two matrix operations needed to complete the algorithm
function ret = matfib(n)
A = [1,1; 1,0];
endfunction
%!assert (matfib(0) == 0);
%!assert (matfib(1) == 1);
%!assert (matfib(2) == 1);
%!assert (matfib(3) == 2);
%!assert (matfib(4) == 3);
%!assert (matfib(5) == 5);
%!assert (matfib(6) == 8);
%!assert (matfib(25) == 75025);
We can expect the second Fibonacci implementation to be faster for two distinct reasons
It's possible to compute matrix powers rather quickly (O(log n)
compared O(n)), and since the fast algorithm is also simple, we
can hope that octave implements it. Since the source to octave is
available, we could actually check this.
Octave is interpreted, so loops are generally slower than matrix operations (which can be done in a single call to an optimized library). This general strategy is called vectorization, and applies in a variety of languages, usually for numerical computations. In particular most PC hardware supports some kind of hardware vector facility.
For the following, you will need either your solutions from last time, or to download timeit.m and tabfib.m. Note that the code here assumes the repetitions come first for timeit, as expected by the version for this lab.
function bench3
timeit(@tabfib, 42, 10000)
timeit(@matfib, 42, 10000)
endfunction