Part 0) Before writing and evaluting our power function we must create a few functions to help us with the process. This includes a trace function, which will be used in our power function. The second function creates a matrix from a list of eigenvalues that we input. The third function uses the eigenvalues of a matrix to calculate the power function we are trying to create through iteration. The third function only works if the inputed matrix is positive semi-definite.
tr <- function(A) {
return(sum(diag(A)))
}
makemyday <- function(eval) {
n <- length(eval)
eval <- sort(eval, decreasing = TRUE)
k <- qr.Q(qr(matrix(rnorm(n^2), n, n)))
return(k %*% diag(eval) %*% t(k))
}
compute.limit <- function(A) {
eig <- eigen(A)
vals <- eig$values
vecs <- eig$vectors
k1 <- vecs[, which.max(vals)] # get the eigenvector for the largest eigenvalue
this <- k1 %*% t(k1)
return(this)
}
Part 1) In order to write a power function that converges we must first work through the process of raising a matrix to the appropriate power.
Bk <- function(B, k) {
# raises a matrix to the k power and divides by the trace of B to the k
# power.
# Args: B: Matix to be raised to a given power. k: the power which B is
# raised to. Returns: Matrix B raised to the kth power, divided by the
# trace of B to the k power.
B1 <- B
i <- 1
while (i < k) {
# Matix multiplies B to itself k-1 times
B1 <- B1 %*% B
i <- i + 1
}
Bk <- B1/tr(B1)
return(Bk)
}
Z <- makemyday(c(4, 54, 6, 3, 4))
Bk(Z, 1)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.27122 0.15824 0.06959 -0.19852 0.18965
## [2,] 0.15824 0.18091 0.04795 -0.13309 0.13302
## [3,] 0.06959 0.04795 0.07690 -0.06878 0.05617
## [4,] -0.19852 -0.13309 -0.06878 0.26086 -0.18389
## [5,] 0.18965 0.13302 0.05617 -0.18389 0.21010
Part 2) Some of the assumptions we will make for this convergence function are as follows:
n.Bkit <- function(B, itmax = 100, eps = 1e-06, verbose = FALSE) {
# finds the matrix the Bk converges to as k goes to invinity.
# Args: B: Matix we wish to converge. itmax: max number of iterations
# allowed. eps: how precise we wish our function to be. verbose: whether
# to print the first five steps of convergence or not. Returns: the matrix
# Bk converges to when k is large.
if (tr(B) == 0) {
stop("The trace of your matrix must be greater than 0.")
}
Bkit <- B/tr(B)
i <- 1
while (i <= itmax) {
if (verbose && i <= 5) {
print(i)
print(Bkit)
}
Bkit2 <- (Bkit %*% B)/tr(Bkit %*% B)
check <- sqrt(sum((Bkit - Bkit2)^2))
Bkit <- Bkit2
i <- i + 1
if (check < eps) {
break
}
}
return(Bkit)
}
We now test our function using a matrix generated by the makemyday function given above (for the purpose of creating matrices that fit our assumptions) and compare it to the compute.limit function given above (for the purpose of testing our code).
Bkit(Z)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.30522 0.21891 0.09743 -0.29051 0.2651
## [2,] 0.21891 0.15700 0.06988 -0.20836 0.1901
## [3,] 0.09743 0.06988 0.03110 -0.09273 0.0846
## [4,] -0.29051 -0.20836 -0.09273 0.27651 -0.2523
## [5,] 0.26505 0.19010 0.08460 -0.25228 0.2302
compute.limit(Z)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.30522 0.21891 0.09743 -0.29051 0.2651
## [2,] 0.21891 0.15700 0.06988 -0.20836 0.1901
## [3,] 0.09743 0.06988 0.03110 -0.09273 0.0846
## [4,] -0.29051 -0.20836 -0.09273 0.27651 -0.2523
## [5,] 0.26505 0.19010 0.08460 -0.25228 0.2302
In the next parts we examine how our function is affected when the assumption we made to start are violated.
Part 3)
Our first example is when our matrix is not positive semi-definite. You will notice across the iterations that [1,1] oscilates as it is calculated. At first it goes from negative to positive and back. As it approaches a convergence point, instead of switching between positive and negative it increases a bit then decreases a bit.
set.seed(11345)
Y <- makemyday(c(-3, 4, 9, -2))
Bkit(Y, verbose = T)
## [1] 1
## [,1] [,2] [,3] [,4]
## [1,] -0.36755 0.04751 -0.0290 0.03758
## [2,] 0.04751 0.01455 0.1273 0.38202
## [3,] -0.02900 0.12730 1.0585 0.14661
## [4,] 0.03758 0.38202 0.1466 0.29451
## [1] 2
## [,1] [,2] [,3] [,4]
## [1,] 0.081223 -0.003552 -0.004935 0.006488
## [2,] -0.003552 0.095777 0.111261 0.080591
## [3,] -0.004935 0.111261 0.674297 0.143073
## [4,] 0.006488 0.080591 0.143073 0.148704
## [1] 3
## [,1] [,2] [,3] [,4]
## [1,] -0.034405 0.006568 -0.00822 0.003346
## [2,] 0.006568 0.053608 0.16472 0.088815
## [3,] -0.008220 0.164716 0.86958 0.212819
## [4,] 0.003346 0.088815 0.21282 0.111221
## [1] 4
## [,1] [,2] [,3] [,4]
## [1,] 0.011684 -0.001146 -0.005593 0.000874
## [2,] -0.001146 0.049107 0.170155 0.062300
## [3,] -0.005593 0.170155 0.853247 0.221705
## [4,] 0.000874 0.062300 0.221705 0.085963
## [1] 5
## [,1] [,2] [,3] [,4]
## [1,] -0.0038421 0.0001483 -0.005806 -0.0005193
## [2,] 0.0001483 0.0426608 0.180858 0.0573590
## [3,] -0.0058057 0.1808576 0.885653 0.2360385
## [4,] -0.0005193 0.0573590 0.236038 0.0755283
## [,1] [,2] [,3] [,4]
## [1,] 3.188e-05 -0.001111 -0.005342 -0.001452
## [2,] -1.111e-03 0.038699 0.186114 0.050607
## [3,] -5.342e-03 0.186114 0.895090 0.243385
## [4,] -1.452e-03 0.050607 0.243385 0.066179
When the largest eigenvalue in modulus is a negative number you get the same result of oscilation however you won't expect it to be switching around a negative number because you have two negatives cancel each other out.
set.seed(11345)
X <- makemyday(c(-100, 2, -4, 30))
Bkit(X, verbose = T)
## [1] 1
## [,1] [,2] [,3] [,4]
## [1,] 1.35289 -0.21034 0.04545 0.03254
## [2,] -0.21034 0.04665 -0.08008 -0.06441
## [3,] 0.04545 -0.08008 -0.37449 -0.09289
## [4,] 0.03254 -0.06441 -0.09289 -0.02506
## [1] 2
## [,1] [,2] [,3] [,4]
## [1,] 0.89139 -0.1424707 0.02767 0.0249385
## [2,] -0.14247 0.0270496 0.01077 -0.0003779
## [3,] 0.02767 0.0107651 0.07470 0.0207705
## [4,] 0.02494 -0.0003779 0.02077 0.0068669
## [1] 3
## [,1] [,2] [,3] [,4]
## [1,] 1.00031 -0.159954 0.03171 0.028269
## [2,] -0.15995 0.024556 -0.01021 -0.005954
## [3,] 0.03171 -0.010211 -0.02384 -0.005860
## [4,] 0.02827 -0.005954 -0.00586 -0.001023
## [1] 4
## [,1] [,2] [,3] [,4]
## [1,] 0.96553 -0.154439 0.030419 0.027241
## [2,] -0.15444 0.025013 -0.003377 -0.003953
## [3,] 0.03042 -0.003377 0.008153 0.002815
## [4,] 0.02724 -0.003953 0.002815 0.001302
## [1] 5
## [,1] [,2] [,3] [,4]
## [1,] 0.97573 -0.156058 0.030797 0.0275445
## [2,] -0.15606 0.024867 -0.005377 -0.0045283
## [3,] 0.03080 -0.005377 -0.001209 0.0002760
## [4,] 0.02754 -0.004528 0.000276 0.0006162
## [,1] [,2] [,3] [,4]
## [1,] 0.97336 -0.155682 0.0307090 0.0274741
## [2,] -0.15568 0.024900 -0.0049117 -0.0043943
## [3,] 0.03071 -0.004912 0.0009687 0.0008668
## [4,] 0.02747 -0.004394 0.0008668 0.0007755
Part 4) In this part we test what happens when there is no unique largest eigenvalue. In particular our estimate will be off because the convergence function will have an extra one in the denominator. \[ \{D_k\}_{ii}=\frac{\left\{\frac{\lambda_i}{\lambda_1}\right\}^k}{1+1+\left\{\frac{\lambda_3}{\lambda_1}\right\}^k+\cdots+\left\{\frac{\lambda_n}{\lambda_1}\right\}^k}. \]
set.seed(11345)
W <- makemyday(c(22, 1, 2, 22))
Bkit(W, verbose = T)
## [1] 1
## [,1] [,2] [,3] [,4]
## [1,] 0.02417 0.019619 -0.013037 0.02320
## [2,] 0.01962 0.180164 0.005986 0.19859
## [3,] -0.01304 0.005986 0.467600 -0.00308
## [4,] 0.02320 0.198588 -0.003080 0.32806
## [1] 2
## [,1] [,2] [,3] [,4]
## [1,] 0.003808 0.019384 -0.014451 0.027489
## [2,] 0.019384 0.164181 0.006833 0.230128
## [3,] -0.014451 0.006833 0.496889 -0.003552
## [4,] 0.027489 0.230128 -0.003552 0.335122
## [1] 3
## [,1] [,2] [,3] [,4]
## [1,] 0.002787 0.019189 -0.014544 0.027906
## [2,] 0.019189 0.162492 0.006904 0.232958
## [3,] -0.014544 0.006904 0.499244 -0.003594
## [4,] 0.027906 0.232958 -0.003594 0.335477
## [1] 4
## [,1] [,2] [,3] [,4]
## [1,] 0.002737 0.01916 -0.014550 0.027944
## [2,] 0.019163 0.16233 0.006910 0.233206
## [3,] -0.014550 0.00691 0.499436 -0.003597
## [4,] 0.027944 0.23321 -0.003597 0.335495
## [1] 5
## [,1] [,2] [,3] [,4]
## [1,] 0.002735 0.01916 -0.014551 0.027947
## [2,] 0.019160 0.16232 0.006910 0.233228
## [3,] -0.014551 0.00691 0.499453 -0.003598
## [4,] 0.027947 0.23323 -0.003598 0.335495
## [,1] [,2] [,3] [,4]
## [1,] 0.002735 0.019160 -0.014551 0.027947
## [2,] 0.019160 0.162315 0.006911 0.233231
## [3,] -0.014551 0.006911 0.499455 -0.003598
## [4,] 0.027947 0.233231 -0.003598 0.335496
compute.limit(W)
## [,1] [,2] [,3] [,4]
## [1,] 0.00497 0.02977 -0.04601 0.04406
## [2,] 0.02977 0.17837 -0.27565 0.26398
## [3,] -0.04601 -0.27565 0.42599 -0.40795
## [4,] 0.04406 0.26398 -0.40795 0.39067
Part 5)
In this example we see when the matrix is of rank p (in our example p = 2 and n = 5 leaving three 0's). What we see is that the matrix will converge instantly to its desired output. However this matrix violates our unique largest eigenvalue assumption and as such compute.limit will give a different result.
V <- makemyday(c(1, 1, 0, 0, 0))
Bkit(V)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.16662 0.03941 -0.190846 -0.13200 0.012165
## [2,] 0.03941 0.42239 -0.061498 0.10357 -0.129314
## [3,] -0.19085 -0.06150 0.219246 0.14586 -0.008699
## [4,] -0.13200 0.10357 0.145858 0.14856 -0.052773
## [5,] 0.01217 -0.12931 -0.008699 -0.05277 0.043193
compute.limit(V)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0.00000 0.000000 0.00000 0.00000
## [2,] 0 0.82613 -0.032718 0.26958 -0.26438
## [3,] 0 -0.03272 0.001296 -0.01068 0.01047
## [4,] 0 0.26958 -0.010676 0.08797 -0.08627
## [5,] 0 -0.26438 0.010471 -0.08627 0.08461