Programming

• Week 1
• Week 2
• Week 3

Week 1: Matlab tutorial and Fibonacci search algorithm

We use the Fibonacci search algorithm to find the approximate minimum of a unimodal nonlinear function.

Step 1: Create the function file (f.m)

1. Choose any nonlinear function that is unimodal in the range [a,b]. You can check to see if your function is unimodal by graphing it over the range as follows:
   
   x = a:1/100:b; %(x has the range [a,b] and we want to plot 100 points between)
   y = x.^2; %(x is automatically created as a matrix, so we must use the operator for matrix arithmetic here, that is, ".^")
   plot(x,y) %(place a semi-colon on the lines you do not wish to print, but leave it off for printed lines)
2. Open a new .M file called f.m - this is just for the function. In Matlab, we describe a function in much the same way as other programming languages. The function can accept input arguments and can output an array of solutions.
3. Enter the function.
   
   %Function file f.m
   function fp = f(p)
   fp = p^2
   
   "function" tells Matlab that you are creating a function called "f" (the name of the file must be the name of the function).
   "fp" is a dummy name for the output of the function.
   "f(p)" tells Matlab you are taking p as an input argument.
   Note that this function is taking one value, p, rather than a matrix. Hence we do not need to use the operator for matrix arithmetic here.

Step 2: Create the Fibonacci number file (fibonacciNumber.m)

This file outputs Fibonacci number n for our use in the Fibonacci search algorithm.

1. (a) Open a new .M file and save it as fibonacciNumber.m.
2. (b) Enter the following text into the file:
   
   % Matlab function m-file: fibonacciNumber.m
   %
   % input: positive integer n
   % output: the n^th number of the Fibonacci sequence
   
   function number = fibonacciNumber(n)
   
   u = sqrt(5);
   number = (1/u)*( ((1+u)/2)^(n+1) - ((1-u)/2)^(n+1) );
   
   where "function" is telling Matlab we are creating another function, where the output is "number" and the input argument is "n".

Step 3: Implement the Fibonacci search algorithm (fibonacciSearch.m)

1. (a) Open a new .M file and save it as fibonacciSearch.m.
2. (b) Go to http://www.ms.unimelb.edu.au/~s620361/programming/fibonacciSearch.m and copy the algorithm to your file.
3. (c) In the Matlab console, type:
   
   fibonacciSearch(@f,a,b,0.01)
   
   where "@f" tells Matlab you want to call function file "f"; you are setting the range to a and b; and you want the tolerance value to be 0.01.

Week 2: Golden section search algorithm

Week 2 (Beginner)

1. Function label. First tell Matlab that you are creating a function with your desired output [xminEstimate,fminEstimate,f_calculations]. The function will be called 'goldenSectionSearch' - remember that this should also be the name of the file. You also need to tell Matlab the input parameters you require, which are (functionToMinimise,a,b,tolerance) where [a,b] is our starting interval and 'tolerance' is the final size of the section in our solution.
   
   function [xminEstimate,fminEstimate,f_calculations] = goldenSectionSearch(functionToMinimise,a,b,tolerance)
   
2. Documentation. Every time you write a piece of code it is advisable to describe the purpose of the code, what it does, the input parameters, the output parameters, and any necessary conditions. This is also an appropriate place to credit the author of the code and the location you retrieved it from.
   
   %GOLDEN Minimize function of one variable using golden section search
   %
   % [xminEstimate,fminEstimate,f_calculations] = goldenSectionSearch(functionToMinimise,a,b,tolerance)
   % computes a local minimum
   % of functionToMinimise. xminEstimate is the computed local minimizer of functionToMinimise and fminEstimate is
   % f(xminEstimate). xminEstimate is computed to an relative accuracy of tolerance.
   %
   % The function must be unimodal on the interval [a,b]:
   % xminEstimate satisfies a < xminEstimate < b. golden is guaranteed to succeed if f
   % is continuous between a and b
   %
   
3. Error checking. It is not necessary, but nonetheless considered to be good practice to perform error checking along the way. We have necessary conditions for this algorithm to work, so here we check that the input parameters satisfy these:
   
   % Check that the correct number of input arguments have been passed to the
   % function (this is not essential, but it helps to detect user error)
   
   if (nargin ~= 4) % in Matlab the operator ~= means "not equal to"
   error('four input arguments are required')
   end
   
   % Check that input parameters have appropriate values
   
   if(b <= a)
   error('b must be strictly greater than a')
   end
   
   if(tolerance <= 0)
   error('tolerance must be strictly positive')
   end
   
4. Algorithm.
   % begin the Golden Search algorithm
   
   gamma = (sqrt(5) - 1)/2;
   
   k = 1; % iteration counter
   
   
   p = b - gamma*(b-a);
   q = a + gamma*(b-a);
   
   fp = feval(functionToMinimise,p);
   fq = feval(functionToMinimise,q);
   f_calculations = 2; % f-calculation counter
   
   while ( (b-a) >= 2*tolerance )
   
   k = k + 1;
   f_calculations = f_calculations + 1;
   
   if (fp <= fq)
   
   b = q;
   q = p;
   fq = fp;
   
   p = b - gamma*(b-a);
   fp = feval(functionToMinimise,p);
   
   else
   
   a = p;
   p = q;
   fp = fq;
   
   q = a + gamma*(b-a);
   fq = feval(functionToMinimise,q);
   
   end
   
   end
   
   % assign the output values of this function
   
   xminEstimate = (a+b)/2; % the midpoint of the final interval
   fminEstimate = feval(functionToMinimise,xminEstimate);
   
5. Printing the solution.
   % Print out the endpoints of the final interval and the output values of
   % the function (can be commented out)
   
   sprintf(' The golden section search algorithm used %d iterations. \n The minimum lies in the interval [%f, %f] \n xminEstimate = midpoint = %f \n fminEstimate %f',k,a,b,xminEstimate,fminEstimate)
   
6. Executing the function.
   Once you have implemented the algorithm, you need to:
   
   1. Create a function file called functionToMinimise.m and save a nonlinear function in it. For example:
      function f = functionToMinimise(x)
      f = x^2 %f and x are dummy variables
      
   2. Go to the main Matlab console and type
      >> goldenSectionSearch(@functionToMinimise, a, b, tolerance) %where you insert the interval [a,b] and final tolerance
      

Week 2 (Advanced)

1. Implement your own code for the Golden Section search algorithm (http://www.ms.unimelb.edu.au/~s620361/material Lecture Friday 7th March).
2. What advantages does the golden section search have over the Fibonacci search?
3. Also implement the false position method and/or Newton's method.

Week 3: A Procedure to Find a Step Size that Satisfies The Armijo-Goldstein and Wolff Conditions

By the time you complete this unit, we want you to be proficient in algorithm implementation - a highly prized skill! This will require some thought and effort on your behalf. Today we will attempt to implement an algorithm from scratch. If you would like to try this independently, turn to page 35 of your lecture notes and start implementing A Procedure to Find a Step Size that Satisfies The Armijo-Goldstein and Wolff Conditions. If you would like some clues on where to start, read the following guidelines (we encourage you to refer to the last two labs for examples if you don't know what something means - but feel free to ask us also!):

1. Open and save. Open a new .M file and save it as 'agw_step.m'. This is where you will create the function that runs the algorithm.
2. Function definition. The first line of your code should be telling Matlab what you are doing.
   1. You are creating a function with an array of outputs (choose them yourself... what would you like to know? the minimum? number of iterations? a reminder of tolerance?)
   2. and the function name 'agw_step' (it MUST match the name of the file)
   3. that accepts certain input parameters (think about what these must be... the function itself? starting point? tolerance? does the algorithm need anything else?)
3. Documentation. The next section of your code should be documentation, which is a commented (%) section that explains the purpose of the code, what the input and outputs mean, and the necessary conditions for the code to work. Sometimes it feels redundant writing documentation, as you may think that no-one but yourself will ever look at your program. The important thing is that if you choose to look at your program in one years time, you can still tell exactly what it does. Think of this work as creating a library of algorithms for your own use later - make sure you have enough information in the documentation for the program to be useful. You should also include your name as the program author, and the date you wrote the program. If you make subsequent improvements to the program you can make notes on these in the documentation.
4. Initialisation.
   1. Initial values. If there are any initial values for the algorithm, you can allocate them to the appropriate parameters here. For example, tlo = 0.
   2. Counter. It is usually interesting to know how many iterations the algorithm performed before it arrived at the solution. You can create a counter variable here, set it to zero, and then add one to it every time an iteration occurs. In the FibonacciSearch algorithm, the counter was called f_calculations.
   3. Calculate function. Before you begin a loop (such as if or while), it is a good idea to calculate any function values that will be needed repeatedly (but do not get updated at each calculation). This way we only have to calculate it once rather than every time we perform an iteration. For example, in this procedure we need to use the value f(0) many times. We only really want to calculate it once, so we set f0 = feval(f,0). This way we can just refer to the parameter f0 whenever we mean f(0).
5. Loops. Algorithms generally involve performing many iterations of a particular calculation, taking into account the previous calculations - this is why we get a computer to do it for us. These steps that are repeated can be thought of as loops. That is, you want the computer to keep going through a list of calculations until some requirement is met. Around and around. It is important, therefore, that these loops eventually come to an end (otherwise it will continue on forever and you will have to force-quit Matlab).
   1. while loops.
      The procedure gives us a condition that is a clue to how to end this looping - it is the last line of the procedure. You can create the loop by telling Matlab:
           while(condition is unmet)
                calculations
           end
      You can collect two conditions together by:
           ((condition1 is unmet)||(condition2 is unmet))
      where || means "or". That is, while one or the other condition is unmet, we want the computer to continue with the calculations.
   2. if conditional statements.
      We often have the condition that if(some condition is met) then (something happens). We sometimes have the "otherwise" condition. In this case we can use else as well. Here is an example:
           if t > 3
                t = t/2
           else
                t = t + 1
           end
      This is saying that if t is greater than 3, then we want to set the new t value as half of old t. Otherwise (if t is less than or equal to 3) we want to set the new t value as old t plus 1. In Matlab you type end to let it know that this particular conditional statement has come to an end.
6. Create function files. Just as you did in the last two labs, create a f.m file and f'.m file (for the derivative).
7. Call the function. In the Matlab main console, type agw_step(@f,@f',0.1,0.6,0.7). Once it works, try different initial t, sigma and mu values.