1.5 Higher Order Functions

One of the things we should demand from a powerful programming language is the ability to build abstractions by assigning names to common patterns and then to work in terms of the names directly. Functions provide this ability.

Also, there are common programming patterns that recur in code, but are used with a number of different functions. These patterns can also be abstracted, by giving them names.

To express certain general patterns as named concepts, we will need to construct functions that can accept other functions as arguments or return functions as values. Functions that manipulate functions are called higher-order functions.

---

1.5.1 Functions as arguments

Suppose we have the following R functions:

1. Sum of natural numbers up to n $$\sum_{k=1}^n k$$

   sum_naturals <- function(n){
       total <- 0
       for(k in 1:n){
           total <- total + k
       }
       return(total)
   }

2. Sum of the cubes of natural numbers up to n.

   $$\sum_{k=1}^n k^3$$

   sum_cubes <- function(n){
       total <- 0
       for(k in 1:n){
           total <- total + k*k*k
       }
       return(total)
   }

3. Sum of the series $\frac{8}{1\times 3} + \frac{8}{5\times 7} + \frac{8}{9 \times 11}\cdots$ which converges to $\pi$

   $$\sum_{k=1}^n \frac{8}{(4k-3)\times(4k-1)}$$

   pi_sum <- function(n){
       total <- 0
       for(k in 1:n){
           total <- total + 8/((4*k-3) * (4*k-1))
       }
       return(total)
   }

These functions share a common underlying pattern. They are the same in most parts, only differing in:

• function name

• the function <func> applied on k inside the loop

<name> <- function(n){

  total <- 0
  for(k in 1:n){
    total <- total + <func>(k)
  }
  
  return(total)
  
} 

If that is the case, we can use a function as an argument input of another function.

Suppose the following summation function, with an additional argument func.

summation <- function(n, func){
    total <- 0
    for(k in 1:n){
        total <- total + func(k)
    }
    return(total)
}

This implements the following summation:

$$\sum_{k=1}^n f(k)$$

Examples

1. Summation of cube

   $$f(k) = k \times k \times k$$

   cube <- function(k){
       k*k*k
   }
   summation(10, cube)

   ## [1] 3025

2. Approximation of $\pi$

   $$f(k) = \frac{8}{(4k-3)\times(4k-1)}$$

   pi_term <- function(k){
       8/((4*k-3)*(4*k-1))
   }

   summation(1000000, pi_term)

   ## [1] 3.141592

   ---

1.5.2 Functions as general methods

With higher-order functions, we begin to see a more powerful kind of abstraction: some functions express general methods of computation, independent of the particular functions they call.

Despite this conceptual extension of what a function means, our environment model of how to evaluate a call expression extends gracefully to the case of higher-order functions, without change.

When a user-defined function is applied to some arguments, the formal parameters are bound to the values of those arguments (which may be functions) in a new local frame.

improve: An iterative improvement algorithm
An iterative improvement algorithm begins with a “guess” of a solution to an equation. It repeatedly applies an update function to improve that guess, and applies a close comparison to check whether the current guess is “close enough” to be considered correct.

improve <- function(update, close, guess=1){
   while(!close(guess)){
    guess <- update(guess)
   }
   return(guess)
}

The improve function is a general expression of repetitive refinement. It doesn’t specify what problem is being solved: those details are left to the update and close functions passed in as arguments.

Example: Golden Ratio

Consider the following example, which implements a general method for iterative improvement and uses it to compute (or estimate) the golden ratio.

The golden ratio $\phi$ , often called “phi”, is a number near 1.6 that appears frequently in nature, art, and architecture.

Among the well-known properties of the golden ratio are that it can be computed by repeatedly summing the inverse of any positive number with 1, and that it is one less than its square.

$$\begin{align} \text{algebraic form}\\ \varphi &=\frac{1+\sqrt{5}}{2}\\ \text{property 1}\\ \varphi &= 1 + \frac{1}{1 + \frac{1}{1 + \ddots}}\\ &= 1+\frac{1}{\varphi} \\ \text{property 2}\\ \varphi &= \varphi^2-1 \end{align}$$

These two properties may be represented by the following functions to be used with improve:

golden_update <- function(guess){
  return(1 + 1/guess)
}

square_close_to_successor <- function(guess){
  return(approx_eq(guess, guess^2 - 1))
}

Above, we introduced the function approx_eq that is used to return TRUE if the two values are approximately equal to each other at some tolerance value.

approx_eq <- function(x, y, tolerance=1e-15){
  return(abs(x - y) < tolerance)
}

Calling improve with the arguments golden_update and square_close_to_successor will compute a finite approximation to the golden ratio.

improve(golden_update, square_close_to_successor)

## [1] 1.618034

Tracing the steps…

1. A local frame for improve is constructed with bindings for update, close, and guess.

2. In the body of improve, the name close is bound to square_close_to_successor, which is called on the initial value of guess.

To visualize without using functions, here is the implementation of the steps used:

# initial guess for phi
guess <-1 

# while property 2 is not satisfied, phi will be updated
repeat{
    # update phi using property 1
    guess <- 1 + 1/guess
    
    # compute phi using property 2
    prop2 <- guess^2 - 1
    
    # exit if guess^2 - 1 and guess are so close to each other
    if (abs(prop2 - guess) < 1e-15){
        break
    }
}
print(guess)

## [1] 1.618034

The value 1.618034 is accurate up to 15 decimal places.

This example illustrates two related big ideas in computer science.

• First, naming and functions allow us to abstract away a vast amount of complexity. While each function definition has been trivial, the computational process set in motion by our evaluation procedure is quite intricate.

• Second, it is only by virtue of the fact that we have an extremely general evaluation procedure for the R language that small components can be composed into complex processes.

Understanding the procedure of interpreting programs allows us to validate and inspect the process we have created.

As always, our new general method improve needs a test to check its correctness.

The golden ratio can provide such a test, because it also has an exact closed-form solution $\varphi=(1+\sqrt{5})/2$, which we can compare to the iterative result of improve

phi <- 1/2 + sqrt(5)/2

approx_phi <- improve(golden_update, square_close_to_successor)

approx_eq(phi, approx_phi)

## [1] TRUE

---

1.5.3 Nested Definitions

Problems we might encounter with what we did earlier:

• Global frame becomes cluttered with names of small functions, which must all be unique.
• We are constrained by particular function signatures: the update argument to improve must take exactly one argument.

Nested function definitions address both of these problems but require us to enrich our environment model.

Example: Square root of a number

Let’s consider a new problem: computing the square root of a number.

Repeated application of the following update converges to the square root of a number $S$.

(Babylonian Method/Heron’s Method):

1. Begin with an arbitrary starting value $x_0$ (close to the actual value of square root of $S$, the better).

2. Let $x_{n+1}$ be the average of $x_n$ and $S/x_n$

3. Repeat step 2 until the desired accuracy is achieved.

It can also be represented as:

$$\begin{align} x_0 &\approx \sqrt{S} \\ x_{n+1} &=\frac{1}{2}\left( x_n +\frac{S}{x_n}\right)\\ \sqrt{S} &=\lim_{n\rightarrow\infty}x_n \end{align}$$

Using the repeat loop, we can use this algorithm to estimate the square root of $S=5$:

S <- 5

# initial guess for sqrt(5)
x <- 2
repeat{
    # update value of square root
    x <- (x + S/x)/2
    # exit if x^2 and S are close
    if(abs(x*x-S) <1e-15){
        break
    }
}

But the goal of this section is not to implement loops, but to use higher order functions, specifically nested definitions.

The repeated application of the following functions converges to the square root of $S$.

average <- function(x, y){
  return((x + y)/2)
}

sqrt_update <- function(x, S){
  return(average(x, S/x))
}

This two-argument update function is incompatible with improve (they both take two arguments, not one), and it provides only a single update, while we really care about taking square roots by repeated updates.

The solution to both of these issues is to place function definitions inside the body of other definitions.

sqrt_cust <- function(S, guess){
    # first function defined inside a function definition
    sqrt_update <- function(x){
    return(average(x, S/x))
    }
    
    # second function defined inside a function definition
    sqrt_close <- function(x){
    return(approx_eq(x*x, S))
    }
    
    
    return(improve(sqrt_update, sqrt_close, guess))
}

sqrt_cust(5, guess = 2)

## [1] 2.236068

Lexical Scope

Locally defined functions also have access to the name bindings in the scope in which they are defined.

This discipline of sharing names among nested definitions is called lexical scoping. Critically, the inner functions have access to the names in the environment where they are defined (not where they are called).

The scoping rules of a language determine how a value is associated with a free variable in a function.

Consider the function below:

func <- function(a, b){
    x + (a*b)
}

Free variables are not formal arguments and are not local variables (assigned inside the function body). In this example, x is a free variable.

R uses lexical scoping or static scoping.

An alternative to lexical scoping is dynamic scoping which is implemented by some languages.

We require two extensions to our environment model to enable lexical scoping.

1. Each user-defined function has a parent environment: the environment in which it was defined.

2. When a user-defined function is called, its local frame extends its parent environment.

To further illustrate, we have another example

y <- 10 
f <- function(x){
    y <- 2
    y^2 + g(x)
}

g <- function(x){
    x*y
}

The y in f()

The y in the definition of f() is not a formal parameter, but defined inside the body. This is a local variable.

The y in g()

On the other hand, in the new function g(), y is still not a formal argument and also not locally defined. This is a free variable.

Its value will be based on the parent environment where it was defined, not on the local frame of g().

Questions:

• What is the value of f(3)?

• What is the value of g(3)?

1.5.4 Functions as return values

Another power of functions is that they can also return another function.

Function composition is a natural method of combination to include in our programming language.

That is, given two functions f(x) and g(x), we might want to define h(x)=f(g(x))

We can define function composition using our existing tools

Notice in the following function compose1, another function is defined instead, which is h.

compose1 <- function(f, g){
  h <- function(x){
    return(f(g(x)))
  }
  return(h)
}

For example, the goal is to compute the square of the input value + 1.

That is: $$(x + 1)^2$$

Let’s apply this function.

square <- function(x){
  return(x*x)
}

successor <- function(x){
  return(x + 1)
}

compose1 <- function(f, g){
  h <- function(x){
    return(f(g(x)))
  }
  return(h)
}

f <- function(x){
  # Never called.
  return(-x)
}

square_successor <- compose1(square, successor)

The square_successor will not have a return value yet, since it became a function as well.

square_successor

## function (x) 
## {
##     return(f(g(x)))
## }
## <environment: 0x0000023533f11888>

Now, applying this new defined function square_successor:

result <- square_successor(12)
result

## [1] 169

---

Practice Exercises

1. Applying a summary function

   In R, simple summary statistics can be obtained using built-in functions such as sum, mean, sd, and median. Write a simple function, apply_function, that takes in two arguments: f, the function to be applied, and data which is a vector that contains the values to which the function will be applied to. To test your function, assign the following values to x.

   x <- c(1, 2, -3, 6, 9)

   Your function apply_function should return the value returned by the function f when applied to the vector data.

2. Applying a summary function on a trimmed vector

   Modify your function apply_function such that it trims (removes) the maximum and minimum values in the vector data. Use the functions min and max and Boolean logic to accomplish this. Name your function apply_trimmed.

   To exclude a value in a vector, you may use indexing, such as the following line of code:

   # remove 2 from the vector x
   x[x != 2]

   ## [1]  1 -3  6  9

3. Applying a function twice

   Write an R function called apply_twice that takes function (f) and a value x as arguments. This function should apply f to x twice in succession and return the result.

   For example, the expected output of apply_twice(square, 2) is 16, which performs $(2^2)^2 = 16$

4. Generating a sequence

   Write a higher-order function in R called generate_sequence that takes a function (f), a starting value (x), and a number of steps (n).

   The function should generate a sequence by repeatedly applying f to the previous result, starting from x, and store all intermediate results in a vector.

   1. Define the higher-order function generate_sequence(f, x, n).

   2. The function should create a vector of length n+1 where

      • The first element is x
      • Each subsequent element is obtained by applying f to the previous element.
      • To create a vector, use the function c().

   3. Return the generated sequence.

   4. Test the function with the different sequences:

   • Adding 3 (arithmetic progression)

     add_three <- function(x) {x + 3}
     generate_sequence(add_three, 0, 5)
     # Expected Output: c(0, 3, 6, 9, 12, 15)

   • Doubling numbers (geometric progression)

     double <- function(x) {x * 2}
     generate_sequence(double, 1, 5)
     # Expected Output: c(1, 2, 4, 8, 16, 32)

5. Getting the first derivative