Midterm 1

require(Matrix)

## Loading required package: Matrix

options(warn=-1)
source("http://www.macalester.edu/~dshuman1/data/365/365Functions.r")

Question 5:

We want to know whether the Jacobian converges faster with greater diagonal dominance. My initial hypothesis is that this does, in fact, happen. I will consider 3 cases, (3) where the dominant nonzero diagonal grows exponentially faster nonzero, dominated, non-diagonal values. (2) The case where everything but the diagonal is zero - but the diagonal grows. (1) The case where only the diagonal grows, but the other elements on the diagonal are nonzero.

This represents an edited version of the Jacobian that only returns steps.

jacobi = function(A,b,m=25,x = rep(0,n),p=2,tol=0.5*10^(-6),history=FALSE) {
  n = length(b)
  if (history) {
    hist = matrix(NA,nrow=length(b),ncol=(m+1))
    hist[,1] = x
  }
  d = diag(A)
  R = A
  R[cbind((1:n),(1:n))] = 0  # allows for r to be sparse  
  steps=0
  for (j in 1:m) {
    x = (b - R %*% x)/d 
    steps = steps+1
    if (history) {hist[,(j+1)] = as.matrix(x)}
    if (vnorm(b-A%*%x,p) <= vnorm(b,p)*tol) break 
  }
  if (history) return(list(x=x,iterations=steps))
  else return(list(x=x,steps=steps))
}

1. Case where the strictly dominant diagonal grows linearly, while the other nonzero entries remain constant. For this, we choose an arbitrary A that meets those criteria.

b = c(1,2,3,4)

TestCase1 = function(a) {
A = rbind(c(a, 2, 3, 4), c(1,a,3,4), c(1,2, a, 4), c(1,2,3,a))
A
}

StoreCase1 = function(a) {
  vector = c(0,100)
  for (i in 1:100) {
    jacob = (jacobi(TestCase1(a),b, m=100))$steps 
    jacob = Matrix(jacob)
    vector[i] = jacob[1,]
    a = a + 1
  }
  vector
  
}

steps = seq(10,109,1)
plot(steps, StoreCase1(10), ylab = "Steps Required", xlab = "Value of Diagonal Entries", main = "Linearly Increasing Diagonal")

Based on the plot, the number of steps required to estimate using the Jacobian decreases logarithmically as the value of the diagonal increases linearly.

2. Case where everything but the diagonal is zero - but the diagonal grows. I will use an all-zero matrix with an increasing diagonal. (Experimenting to see if any funny behavior with zeroes on non-diagonal). Note: Sparse matrix is not used in order to work with my modified Jacobian function.

b = c(4,3,2,1)
TestCase2 = function(a) {
  A = rbind(c(a,0,0,0), c(0,a,0,0), c(0,0,a,0), c(0,0,0,a))
  A
  }

StoreCase2 = function(a) {
  vector = c(0,100)
  for (i in 1:100) {
    jacob = (jacobi(TestCase2(a),b, m=100))$steps 
    jacob = Matrix(jacob)
    vector[i] = jacob[1,]
    a = a + 1
  }
  vector
  
}

steps = seq(10,109,1)
plot(steps, StoreCase2(10), ylab = "Steps Required", xlab = "Value of Diagonal Entries", main = "Zero-Matrix with Increasing Diagonal")

In the case of a zero-matrix with only entries on the diagonal, the matrix cannot converge faster because it is already essentially ‘converged’. I’m showing this case to illustrate the possibility that a function could converge at the onset, meaning more steps cannot improve it - no matter how strictly diagonally dominant it continues to get.

3. Case where the entire matrix is nonzero and the diagonal is strictly dominant. However, the diagonal increases exponentially while everything else increases less-exponentially. In essence, both the diagonal and the non-diagonal grow, but the diagonal grows faster.

b = c(17,12,5,3)
TestCase3 = function(a) {
 A = rbind(c(10^(a+1), 2^a, 3^a, 4^a), c(1^a,10^(a+1),3^a,4^a), c(1^a,2^a, 10^(a+1), 4^a), c(1^a,2^a,3^a,10^(a+1)))
 #A = rbind(c(a, 2, 3, 4), c(1,a,3,4), c(1,2, a, 4), c(1,2,3,a)) #--- this is used to test that my TestCase is still working properly. It is!
  A
  }

StoreCase3 = function(a) {
  vector = c(0,50)
  for (i in 1:50) {
    jacob = (jacobi(TestCase3(a),b, m=100))$steps 
    jacob = Matrix(jacob)
    vector[i] = jacob[1,]
    a = a + 1
  }
  vector
  
}


plot(StoreCase3(10), ylab = "Steps Required", xlab = "Power of Diagonal (other entries, subtract 1)", main = "Zero-Matrix with Increasing Diagonal")

In this situation, it decreases from 2 steps to 1 step once 10^3 becomes the diagonal value.

Given all of our results, I confirm the conclusion that greater strict diagonal dominance causes the Jacobian to converge faster. The only caveat is in the wording - there is the possibility that it converges after the first iteration if the diagonal is large enough - or - if the non-diagonal entries are 0. However, though it cannot possible converge faster, per say, it still converges nonetheless. So in closing, it always converges faster given greater strict dominance when it does not already show convergence in the first iteration.

Question 6

Part (a)

W = spMatrix(16,16, i = c(2:16), j=c(1:15), x=(rep(1,15))) + spMatrix(16,16, i=c(1:15), j=c(2:16), x=(rep(1,15)))
W[1,9] = .5
W[1,8] = 1
W[8,1] = 1
W[9,1] = .5
W[9,8] = 0
W[8,9] = 0
W[8,10] = (5/12)
W[10,8] = (5/12)
W[16,9] = 1
W[9,16] = 1

image(W)

Part (b)

D = diag(rowSums(W))
L = D-W
eigen = eigen(L, only.values = TRUE)
eigen$values

##  [1]  4.451e+00  4.000e+00  3.671e+00  3.516e+00  3.414e+00  3.414e+00
##  [7]  2.313e+00  2.162e+00  2.000e+00  2.000e+00  9.709e-01  6.342e-01
## [13]  5.858e-01  5.858e-01  1.159e-01 -2.477e-16

Part (c)

J = spMatrix(16,16, i = c(2:16), j=c(1:15), x=(rep(1,15))) + spMatrix(16,16, i=c(1:15), j=c(2:16), x=(rep(1,15)))
#This declares the sparse matrix of 1s on the 1st submajor and subminor diagonals.

#FiedlerValue takes a value of a and adds the corresponding a (and other extra entries) to the sparse matrix. It returns the smallest, positive nonzero eigenvalue.
fiedler.value = function(a){
J[1,9] = a
J[1,8] = 1
J[8,1] = 1
J[9,1] = a
J[9,8] = 0
J[8,9] = 0
J[8,10] = (4-(3*a))/6
J[10,8] = (4-(3*a))/6
J[16,9] = 1
J[9,16] = 1
D = diag(rowSums(J))
L = D-J
eigen =eigen(L, only.values = TRUE)
eigenMat = Matrix(eigen$values) #Made 1-column matrix of values. Was the only way I could pick out the second-to-last value without receiving an error. 

n = length(eigenMat)
return (eigenMat[(n-1),]) #because R orders in descending order by default, the second-to-last eigenvalue will be the Fiedler Value (last is zero). 
}


fiedler.plot = function(b){
  vector = c(0,101)
  for (i in 1:101) {
    vector[i] = fiedler.value(b)
    b = b + (4/300)
  }
 vector #returns vector of eigenvalues for all a. 
}
b = seq(0,(4/3), (4/300))
plot(b, fiedler.plot(0), xlab = "Weight (a)", ylab = "Fiedler Value", main = "Fiedler Value vs. Weight (a)")

Part (d): The discrete derivative.

options(digits = 3) 

Tannenbaum = function(c,delta=(4/300)){ #This is the derivative function, which calls the fiedler value function. Named Tannenbaum because it was like Xmas when I got it to work. 
      ((fiedler.value(c+delta) - (fiedler.value(c))))/(delta)}
                                 
g = function(c){ #This function records a list of the derivatives. 
  Dvector = c(0,101)
 
  for (i in 1:101){
    Dvector[i] = Tannenbaum(c) #iterates 101 times through my Derivative function. 
    #stores in vector Dvector
      c = c + (4/300)
      
      }

Dvector #returns the vector of derivative. 

}
c = seq(0,(4/3), (4/300))
plot(c,g(0), xlab = "Weight (a)", ylab = "Derivative", main = "Fiedler Derivative vs. Weight (a)")

abline(h=0)  

Part(e): Finding the root through bisection.

interval = c(0.7,0.9)
bisect(g,interval) #using the professor's bisect function. 

## $root
## [1] 0.796
## 
## $approximations
##  [1] 0.800 0.750 0.775 0.788 0.794 0.797 0.795 0.796 0.796 0.796 0.796
## [12] 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796
## [23] 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796
## 
## $left.brackets
##  [1] 0.700 0.700 0.750 0.775 0.788 0.794 0.794 0.795 0.796 0.796 0.796
## [12] 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796
## [23] 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796
## 
## $right.brackets
##  [1] 0.900 0.800 0.800 0.800 0.800 0.800 0.797 0.797 0.797 0.796 0.796
## [12] 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796
## [23] 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796 0.796
## 
## $fleft
##  [1] 6.79e-03 6.79e-03 3.20e-03 1.46e-03 6.05e-04 1.78e-04 1.78e-04
##  [8] 7.11e-05 1.78e-05 1.78e-05 4.46e-06 4.46e-06 1.13e-06 1.13e-06
## [15] 2.95e-07 2.95e-07 8.70e-08 8.70e-08 3.50e-08 8.99e-09 8.99e-09
## [22] 2.49e-09 2.49e-09 8.67e-10 5.46e-11 5.46e-11 5.46e-11 5.46e-11
## [29] 3.90e-12 3.90e-12 3.90e-12 3.90e-12
## 
## $fright
##  [1] -6.93e-03 -2.48e-04 -2.48e-04 -2.48e-04 -2.48e-04 -2.48e-04 -3.55e-05
##  [8] -3.55e-05 -3.55e-05 -8.86e-06 -8.86e-06 -2.20e-06 -2.20e-06 -5.37e-07
## [15] -5.37e-07 -1.21e-07 -1.21e-07 -1.70e-08 -1.70e-08 -1.70e-08 -4.01e-09
## [22] -4.01e-09 -7.58e-10 -7.58e-10 -7.58e-10 -3.52e-10 -1.48e-10 -4.69e-11
## [29] -4.69e-11 -2.14e-11 -9.05e-12 -2.39e-12

The bisection method determined that the derivative is zero (i.e. fiedler value is maximized) when the weight (a) is ~0.796.