1 function [x, k] = steepest_descent(g, ...

2 x0, ...

3 step_size, ...

4 tolerance, ...

5 max_iterations)

6 %

7 % An implementation of the steepest-descent algorithm, with the

8 % search direction d_{k} = -\nabla f(x_{k}).

9 %

10 % We should terminate when either,

11 %

12 % a) The infinity-norm of the gradient at x_{k} is greater than

13 % ``tolerance``.

14 %

15 % b) We have performed ``max_iterations`` iterations.

16 %

17 % INPUT:

18 %

19 % * ``g`` - the gradient of the function ``f`` we're minimizing.

20 %

21 % * ``x0`` - the starting point for the search.

22 %

23 % * ``step_size`` - a function of x which returns the optimal step

24 % size alpha_{k}. This will generally require more information

25 % than just x; for example it might need the function ``f`` or the

26 % gradient ``g``. However, our caller has this information so it

27 % should be incorporated into the step size algorithm that we

28 % receive.

29 %

30 % * ``tolerance`` - the stopping tolerance. We stop when the norm

31 % of the gradient falls below this value.

32 %

33 % * ``max_iterations`` - a safety net; we return x_{k}

34 % unconditionally if we exceed this number of iterations.

35 %

36 % OUTPUT:

37 %

38 % * ``x`` - the solution at the final value of x_{k}.

39 %

40 % * ``k`` - the value of k when we stop; i.e. the number of

41 % iterations.

42 %

43 % NOTES:

44 %

45 % A specialized implementation for solving e.g. Qx=b can avoid one

46 % matrix-vector multiplication below.

47 %

49 % The initial gradient at x_{0} is not supplied, so we compute it

50 % here and begin the loop at k=0.

51 k = 0;

52 xk = x0;

53 gk = g(xk);

55 while (k <= max_iterations && norm(gk, 'inf') > tolerance)

56 % Loop until either of our stopping conditions are met. This

57 % should also work when x0 is a satisfactory point.

59 dk = -gk;

60 alpha_k = step_size(xk);

61 x_next = xk + (alpha_k * dk);

63 % We potentially just performed one more iteration than necessary

64 % in order to simplify the loop. Note that due to the structure of

65 % our loop, we will have k > max_iterations when we fail to

66 % converge.

67 k = k + 1;

68 xk = x_next;

69 gk = g(x_next);

70 end

72 x = xk;

73 end