Instructions

Exercise 1: Matrix Multiplication

The usual (not element-wise) matrix multiplication \(C = A \cdot B\) is only defined when the number of rows from \(B\) match the number of columns from \(A\). For our purposes, let us say \(A\) has \(m\) rows and \(n\) columns and \(B\) has \(n\) rows and \(p\) columns. We can denote this in a simplified manner by writing \(A\in\mathbb{R}^{m\times n}\) and \(B\in\mathbb{R}^{n\times p}\). The outcome of matrix multiplication is another matrix \(C \in\mathbb{R}^{m\times p}\) which combines the entries of \(A\) and \(B\).

Matrix multiplication can be viewed in several ways:

  1. The entry \((i,k)\) of \(C\), denoted by \(C_{ik}\) (\(i\)-th row, \(k\)-th column) is given by the scalar sum of products: \[C_{ik} = \sum_{j = 1}^{n} A_{ij} \cdot B_{jk}\]

  2. The \(k\)-th column of the product matrix \(C\) can be thought of as the sum of \(A\)’s columns weighted by the columns of \(B\) \[C_{,k} = \sum_{j=1}^{n} B_{jk} \cdot A_{,j},\] where the notation \(C_{,k}\) is used to indicate column \(k\) to match R’s index notation C[,k]. That is, \(C_{,k}\) is a linear combination of \(A\)’s columns and the weights \(B_{jk}\) tell us how much to “mix” for each column. Note that here we have a scalar-vector product.

1(A) Write a function that implements Formula 1 as given above. Call it matmul1(A, B). Your function should take two arguments, A and B, and return their matrix product C. You may assume that A and B are compatible. Use only arithmetic operators and loops, and avoid functions like sum() or apply(). This version requires the most for loops.

matmul1= function(A,B) {
  C= matrix(0, nrow= nrow(A), ncol= ncol(B))
  for(i in 1: (nrow(A))){
    for(k in 1:ncol(B)){
      count=0
       for(j in 1:ncol(A)){
         count= count + A[i,j]*B[j,k] 
          
       } 
      C[i,k]= count 
    }
  }
  return(C)
}

A= matrix(c(1,2,3,4,5,6,7,8), nrow = 2, ncol = 4)
B= matrix(c(3,6,9,12,15,18,21,24), nrow = 4, ncol = 2)
matmul1(A,B)
##      [,1] [,2]
## [1,]  150  342
## [2,]  180  420
A%*%B
##      [,1] [,2]
## [1,]  150  342
## [2,]  180  420

1(B) Write a function that implements Formula 2 as given above. Call it matmul2(A, B). Your function should take two arguments, A and B, and return their matrix product C. You may assume that A and B are compatible. Use only arithmetic operators and loops, and avoid functions like sum() or apply(). This version requires fewer for loops.

matmul2= function(A,B) { 
  C= matrix(0, nrow(A), ncol(B))
  for(k in 1: ncol(B)){
    count=0
    for(j in 1: ncol(A)) {
      
      count= count+ B[j,k]*A[,j]
    }
    C[,k]= count
  }
  return(C)
}
A= matrix(c(1,2,3,4,5,6,7,8), nrow = 2, ncol = 4)
B= matrix(c(3,6,9,12,15,18,21,24), nrow = 4, ncol = 2)
matmul2(A,B)
##      [,1] [,2]
## [1,]  150  342
## [2,]  180  420
A%*%B
##      [,1] [,2]
## [1,]  150  342
## [2,]  180  420

1(C) Now we are going to compare the computation time of our functions matmul1 and matmul2 to R’s implementation of matrix multiplication A %*% B. To compare times, we measure time with the system.time() function. To use system.time(), simply pass an R expression as an argument to the function; for example system.time(A %*% B) will report timing results for running A %*% B.

Use the system.time() function to compare the computation time for multiplying matrices with matmul1(), matmul2(), and A %*% B and save the results in a named list with names matmul1, matmul2, and baseR, respectively. Display the results as a data.frame object.

Look at the elapsed measurements. Which method was the fastest for you? Which one was the slowest?

# DO NOT MODIFY
set.seed(1000)
A <- matrix(rnorm(1000 * 500), nrow = 1000, ncol = 500)
B <- matrix(rnorm(500 * 800), nrow = 500, ncol = 800)
# END DO NOT MODIFY

## Your code here:
##
##
as.data.frame(list(matmul1=system.time(matmul1(A,B)),matmul2= system.time(matmul2(A,B)), baseR= system.time(A%*%B)))
##            matmul1 matmul2 baseR
## user.self    20.99    1.43  0.17
## sys.self      0.05    0.28  0.00
## elapsed      21.21    1.75  0.17
## user.child      NA      NA    NA
## sys.child       NA      NA    NA

Answer here: baseR is fastest, while matmul1 is the slowest computation speed.

Exercise 2: Using while loops in iterative algorithms

The ancient Babylonians knew how to approximate the hypotenuse of a triangle, which involves computing the square root of a non-negative number. The Greek mathematician Heron would later define an explicit algorithm for computing the square root of an arbitrary non-nengative number \(c\), which is often attributed to the Babylonians for their geometric skills. The method turns out to be a special case of Newton’s method, a standard algorithm for minimizing/maximizing functions (e.g. \(\min f(x)\)). The so-called Babylonian method estimates square roots using the iterative recipe:

\[ x_{n+1} = \frac{1}{2} \left[x_{n} + \frac{c}{x_{n}}\right].\] This is an iterative algorithm in which the next iterate (left-hand side), \(x_{n+1}\), is defined by some function of the previous iterates (just \(x_{n}\) in this case) plus some additional problem data (such as the number, \(c\)).

Source: Babylonian Method (Wikipedia)

Credit: Algorithms from THE BOOK (Chapter 1, Section 3)

2(A) Write a function named my.sqrt() that accepts four arguments, given in the following order:

my.sqrt= function(c,x0, maxiter=10, epsilon= .001){
  count=0
  xn= x0
  eps= epsilon-.5
  while(count<= maxiter & eps<epsilon) {
  
    xn1= .5*(xn+(c/xn))
    eps= xn1^(1/2)-c
    xn= xn1
    count=count+ 1
  }
   return(xn)
}
my.sqrt(c = 30,maxiter=500, x0= 3)
## [1] 5.477226
sqrt(30)
## [1] 5.477226
  1. a number c, whose square root we wish to compute,
  2. an initial guess x0, which must be a positive number,
  3. a maximum number of iterations (or steps) maxiter with default value 10, and
  4. a control parameter epsilon with default value 0.001.

This function should implement the Babylonian method, using a single while loop. Use the following outline to guide you.

Outline

Before the while loop

  • Use x0 to initialize your current estimate of the square root.

  • Decide how you’re going to exit the while loop. What variables might you need to define?

  • Provide brief comments for any variables that you define here, as needed.

In the while loop

  • Implement the algorithm map using the formula above.

  • Check the convergence condition \(|x_{n+1}^2 - c| < \epsilon\). Here \(\epsilon\) is the Greek letter corresponding to epsilon. There are many ways to do this check, but make sure that you do not create an infinite loop by accident, or exit early. If the convergence condition is

    • TRUE, you should exit the while loop.

    • FALSE, you should continue the while loop.

  • Check that you have not taken more than the maximum number of iterations maxiter:

    • if TRUE, you should exit the while loop.

    • if FALSE, you should continue the while loop.

  • Update any additional variables that may need it.

After the while loop

  • Before finishing the function, print the number of iterations performed by the algorithm. For example, if the algorithm iterated 5 times, print Iteration Total = 5.

  • The function must return our estimate of \(\sqrt{c}\).

2(B) [Convergence Check] Call your function with c = pi^2, x0 = 1 and maxiter = 100. The algorithm should finish in 5 steps.

my.sqrt= function(c,x0, maxiter=100, epsilon= .001){
  count=0
  xn= x0
  eps= epsilon-.5
  while(count<= maxiter & eps<epsilon) {
  
    xn1= .5*(xn+(c/xn))
    eps= xn1^(1/2)-c
    xn= xn1
    count=count+ 1
  }
   return(xn)
}
my.sqrt(c = pi^2,maxiter=100, x0= 1)
## [1] 3.141593

2(C) [Iterations Check] Call your function with c = pi^2, x0 = 1, and maxiter = 3. The algorithm should only take 3 steps.

my.sqrt(c = pi^2,maxiter=3, x0= 1)
## [1] 3.141757

2(D) Create a new function called my.sqrt.iterates(), which has the same inputs as your function my.sqrt(), but, instead of returning only the estimator of \(\sqrt{c}\), it returns a vector with all of the iterated values. For example, if for a given set of inputs, the function performs 5 iterations, we must return a vector x = c(x1, x2, x3, x4, x5).

2(E) Run iterates <- my.sqrt.iterates(c = 400, x = 1000, maxit = 1000, epsilon = 0.001). Then, create a convergence plot. The plot must include points corresponding to each iteration step of our estimates of \(\sqrt{c}\), , a curve connecting the points, and a horizontal red line at the true square root. Give appropriate title, and labels to the plot:

  1. Label for x-axis must be "Iteration Index".
  2. Label for y-axis must be "Iteration Value"
  3. Title for plot must be "Babylonian Method: Convergence"

Add a legend on the top-right of the plot describing the lines and points of the iteration, and the horizontal red line.

Exercise 3: Two Ways to Calculate the Fibonacci Sequence

The Fibonacci sequence is a classic example of a sequence of numbers \(\{F_0,F_1,F_2,F_3,...\}\) defined through a recursive formula. The sequence of numbers starts with \(F_0 = 0\), \(F_1 = 1\). The following numbers in the sequence are defined as: \[F_{t+2} = F_{t+1} + F_{t}.\] For example, \(F_2 = 2\), \(F_3 = 3\), \(F_4 = 5\), \(F_5 = 8\) and \(F_6 = 13\). In this exercise, you will create two functions that calculate the Fibonacci sequence.

3(A) Create a function my_fibonacci1(t) that, for any positive type integer value t, it calculates the Fibonacci value \(F_t\). To prevent potential errors, use conditionals to verify that the input value t is both an integer type, and positive. If the value is not integer or not positive, print "Error: input does not conform to required format" and return -1.

If the input t is appropriate, move forward with the calculation of the Fibonacci \(F_t\). For this, start with \(F_1 = 1\), \(F_2 = 1\) and use a for() loop to calculate the values \(F_3,F_4,...,F_{t-1}, a_t\) in sequence. Then, return the value \(F_t\).

To verify your function, save in a vector myfibo1 of length 25 the Fibonacci numbers calculated with the function my_fibonacci1() for the values t = 1L,2L,3L,..., 25L. Print the vector myfibo1.

3(B) Notice that, in the previous function, it was necessary to calculate all Fibonacci numbers from 1 to \(t\) to recover the value \(F_t\). This means that the time it requires to calculate \(t\) is increasing with respect to \(t\). Mathematicians in the XIX century realized that it was possible to provide a non-recursive formula for the Fibonacci sequence, that allowed them to recover the value of \(F_t\) without the need to calculate all of \(F_1,F_2,...,F_{t-1}\). The formula is given as follows: \[ F_t = \frac{1}{2} \left( \varphi^t - \big(-\frac{1}{\varphi}\big)^t \right),\] where the value \(\varphi\) is the golden ratio \(\varphi := \frac{1 + \sqrt{5}}{2}\).

Create a function named my_fibonacci2(t) that uses this formula to calculate the Fibonacci numbers. First, just as in part 3(a), do a check of whether the input is of type integer and positive. If the conditions are met, apply the formula above to calculate the \(t\)-the Fibonacci number, and return it.

To verify your function, save in a vector myfibo2 of length 25 the Fibonacci numbers calculated with the function my_fibonacci2() for the values t = 1L,2L,3L,..., 25L. Print the vector myfibo2.