UNB/ CS/ David Bremner/ teaching/ cs2613/ labs/ Lab 19

# Before the lab

## Linear Algebra

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

# Background

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.

# Running Octave

Time
5 minutes
Activity
Demo / discussion
• 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)
``````

# Recursive Fibonacci

Time
10 minutes
Activity
Demo/Group programming.

In L11 we discovered that caching (also known as memoization) could make some recursive functions much faster. We will (re)consider the same example in Octave. Here is a JavaScript recursive function for Fibonacci (slightly modified from L11).

```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/L19/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.

# Table based Fibonacci

Time
25 minutes
Activity
Programming puzzle

We saw in Lab 11 saving previously computing results can give big speedups. The approach of Lab 11 still incurs the overhead of recursive function calls, which in some languages is quite expensive. A more problem specific approach (sometimes called dynamic programming) is to fill in values in a table.

Save the following in `~/cs2613/labs/L19/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);
```

• What are two important differences about array access in Octave compared to Python and JavaScript?

• What is a difference with racket vectors not related to brackets.

# Performance comparison

Time
10 minutes
Activity
Demo / discussion

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?

# Using a constant amount of storage.

Time
25 minutes
Activity
Programming puzzle

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);
```
• This algorithm is quite tricky. To understand how to fill in the blanks
• what happens for `n==0` (first blank)
• what happens for `n==1` (second blank)
• in loop repetition `i`, b is assigned to `(i+2)`nd Fibonacci number.
• Check the change in performance with the following. Although not as spectacular as repacing recursion (in Octave), there should be a noticable speedup.
```function bench2
timeit(@tabfib, 42, 10000)
timeit(@abfib, 42, 10000)
endfunction
```

# Fibonaccci as matrix multiplication

Time
25 minutes
Activity
Implement formula in Octave

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) ]
```

Since matrix exponentiation is built-in to octave, this is particularly easy to implement in octave

Save the following as `~/cs2613/labs/L19/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.

```function bench3
timeit(@abfib, 42, 10000)
timeit(@matfib, 42, 10000)
endfunction
```