Write a SAS Macro and R function to perform the following task.
For a given distribution with two parameters (it can be normal, uniform, or gamma), and an integer \(n\):
Generate an \(n \times n\) symmetric random matrix (say \(A\)). In this matrix, the elements are random variable from the given distribution. Notice that the matrix must be symmetric.
Compute \(A^{10}\). You should think about the best way to compute this. Try to save the computation time. Hint: use eigenvalue decomposition.
I used some functions from the base package, and made a successful attempt to not use additional libraries such as Matrix. While I could have manually built a function to compute the \(VDV^{-1}\) decomposition of a matrix, I felt my use of svd was appropriate. I initially wrote code to format a symmetric matrix, but the code used and cited below was very elegant and significantly shorter than my approach.
Surprisingly, the base package of R does not have a method to compute the power of a matrix. This means to verify my work would require a bit of work from the user. I do not recommend checking my results with other libraries, especially since other approaches use \(VDV^{-1}\) decomposition to compute powers of a matrix (as opposed to traditional matrix multiplication). I have attempted to verify the accuracy of my work in the final section of this page.
set.seed(30)First, let us create a function n_by_n_sym_rand_matrix that fits the description.
n_by_n_sym_rand_matrix = function(dist, n, para_1, para_2){
# This function takes the name of a distribution as a string
# to compute n-squared random variables given two parameters
# acceptable inputs are "normal", "uniform", and "gamma"
if (dist == "normal") {
values = rnorm(n = n**2, mean = para_1, sd = para_2)}
if (dist == "uniform") {
values = runif(n = n**2, min = para_1, max = para_2)}
if (dist == "gamma") {
values = rgamma(n = n**2, shape = para_1, rate = para_2)}
A = matrix(data = values ,nrow = n, ncol = n)
A[lower.tri(A)] = t(A)[lower.tri(A)]
# The above line of code was found at the following link:
# http://stackoverflow.com/a/21825946
# Credit to StackOverflow user: user3318600
return(A)
}Now, we create a second function fast_mat_exp which will raise the matrix to an inputted power.
fast_mat_exp = function(mat, exp) {
# This function takes a matrix and raises it to the exp power
# This is done using the svd function to obtain the matrix's
# VDV^-1 decomposition.
decomp = svd(mat)
D = diag(decomp$d)
result = decomp$u %*% D**exp %*% solve(decomp$u)
return(result)
}Test your SAS Macro and R function with the following set up:
A_1 = n_by_n_sym_rand_matrix(dist = "normal", n=10, para_1=1, para_2=4**.5)
round(A_1,2) # preview the rounded A_1 matrix to two decimal places## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] -1.58 -1.05 0.68 -2.45 -0.99 -0.70 2.95 3.23 2.70 2.13
## [2,] -1.05 -2.64 2.51 2.23 3.34 1.25 1.80 -0.48 -0.13 0.56
## [3,] 0.68 2.51 -0.82 2.45 0.04 -0.24 1.70 5.02 0.70 0.06
## [4,] -2.45 2.23 2.45 0.92 -2.34 -2.43 0.42 0.95 0.14 2.15
## [5,] -0.99 3.34 0.04 -2.34 3.27 2.43 4.15 0.83 5.44 -1.59
## [6,] -0.70 1.25 -0.24 -2.43 2.43 -0.83 2.69 2.00 -3.07 1.52
## [7,] 2.95 1.80 1.70 0.42 4.15 2.69 -0.84 -0.45 -1.13 1.20
## [8,] 3.23 -0.48 5.02 0.95 0.83 2.00 -0.45 3.21 4.14 1.54
## [9,] 2.70 -0.13 0.70 0.14 5.44 -3.07 -1.13 4.14 -0.46 1.42
## [10,] 2.13 0.56 0.06 2.15 -1.59 1.52 1.20 1.54 1.42 -1.77A_1_to_10 = fast_mat_exp(mat=A_1, exp=10) #construct A^10
# preview the first 5 rows and 5 columns of A^10 with scientific notation
format(A_1_to_10[1:5,1:5], scientific=TRUE, digits = 3) ## [,1] [,2] [,3] [,4] [,5]
## [1,] " 5.35e+09" " 2.28e+09" " 5.61e+09" " 8.47e+08" " 9.15e+09"
## [2,] " 2.28e+09" " 2.38e+09" " 3.37e+09" "-3.43e+08" " 5.14e+09"
## [3,] " 5.61e+09" " 3.37e+09" " 6.89e+09" " 4.67e+08" " 1.03e+10"
## [4,] " 8.47e+08" "-3.43e+08" " 4.67e+08" " 6.35e+08" " 4.35e+08"
## [5,] " 9.15e+09" " 5.14e+09" " 1.03e+10" " 4.35e+08" " 1.79e+10"Now, let us test the accuracy of our fast_mat_exp function. We will do this by manually computing \(A_1^{10}\) and then comparing the rounded entries to two decimal places. This is necessary to deal with errors that might rise from floating point.
A_1_to_5_manual = A_1 %*% A_1 %*% A_1 %*% A_1 %*% A_1
A_1_to_10_manual = A_1_to_5_manual %*% A_1_to_5_manual
round(A_1_to_10_manual,2) == round(A_1_to_10,2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [5,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [6,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [7,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [8,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [9,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [10,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEA_2 = n_by_n_sym_rand_matrix(dist = "uniform", n=15, para_1=0, para_2=3)
# preview the first 10 rows and 10 columns rounded A_2 matrix to two decimal places
round(A_2[1:10,1:10],2) ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 2.09 2.75 0.33 0.60 1.22 2.17 0.20 2.13 2.61 1.42
## [2,] 2.75 1.79 0.98 2.09 0.37 0.33 1.62 0.85 0.71 0.03
## [3,] 0.33 0.98 1.97 2.77 1.35 1.48 0.34 0.50 1.32 1.35
## [4,] 0.60 2.09 2.77 0.09 2.23 0.64 0.55 1.86 1.68 0.54
## [5,] 1.22 0.37 1.35 2.23 0.42 1.30 0.12 0.41 2.31 0.10
## [6,] 2.17 0.33 1.48 0.64 1.30 0.06 2.66 0.16 2.67 1.78
## [7,] 0.20 1.62 0.34 0.55 0.12 2.66 0.31 2.91 2.17 1.67
## [8,] 2.13 0.85 0.50 1.86 0.41 0.16 2.91 0.14 0.56 1.06
## [9,] 2.61 0.71 1.32 1.68 2.31 2.67 2.17 0.56 0.41 1.09
## [10,] 1.42 0.03 1.35 0.54 0.10 1.78 1.67 1.06 1.09 0.09A_2_to_10 = fast_mat_exp(mat=A_2, exp=10) #construct A^10
# preview the first 5 rows and 5 columns of A^10 with scientific notation
format(A_2_to_10[1:5,1:5], scientific=TRUE, digits = 3) ## [,1] [,2] [,3] [,4] [,5]
## [1,] "9.23e+11" "8.59e+11" "1.33e+12" "1.29e+12" "1.13e+12"
## [2,] "8.59e+11" "7.99e+11" "1.23e+12" "1.20e+12" "1.05e+12"
## [3,] "1.33e+12" "1.23e+12" "1.90e+12" "1.86e+12" "1.62e+12"
## [4,] "1.29e+12" "1.20e+12" "1.86e+12" "1.82e+12" "1.58e+12"
## [5,] "1.13e+12" "1.05e+12" "1.62e+12" "1.58e+12" "1.38e+12"Again, we will test the accuracy of our computation. We will do this by manually computing the \(A_2^{10}\) and then comparing the rounded entries to one decimal place. Since these matrices are very large, let us count the entries that are equal, and present it as a present it as a table.
A_2_to_5_manual = A_2 %*% A_2 %*% A_2 %*% A_2 %*% A_2
A_2_to_10_manual = A_2_to_5_manual %*% A_2_to_5_manual
results_2 = (round(A_2_to_10_manual,1) == round(A_2_to_10,1))
table(results_2)## results_2
## FALSE TRUE
## 7 218We see even with rounding to one decimal place, we still have some entries that are off by a bit.
A_3 = n_by_n_sym_rand_matrix(dist = "gamma", n=20, para_1=1, para_2=1)
# preview the first 10 rows and 10 columns rounded A_3 matrix to two decimal places
round(A_3[1:10,1:10],2) ## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 2.95 1.14 1.46 2.85 0.07 0.49 2.64 0.72 0.02 2.57
## [2,] 1.14 1.41 1.28 1.46 0.06 2.91 4.75 0.30 1.42 0.05
## [3,] 1.46 1.28 1.76 0.45 0.81 0.10 2.46 0.20 3.59 1.30
## [4,] 2.85 1.46 0.45 1.48 1.43 0.23 0.00 0.58 0.09 2.08
## [5,] 0.07 0.06 0.81 1.43 0.91 2.79 0.12 0.89 1.39 0.60
## [6,] 0.49 2.91 0.10 0.23 2.79 4.70 1.53 5.26 0.23 0.08
## [7,] 2.64 4.75 2.46 0.00 0.12 1.53 0.35 0.39 0.93 0.06
## [8,] 0.72 0.30 0.20 0.58 0.89 5.26 0.39 0.47 0.06 1.02
## [9,] 0.02 1.42 3.59 0.09 1.39 0.23 0.93 0.06 0.75 3.27
## [10,] 2.57 0.05 1.30 2.08 0.60 0.08 0.06 1.02 3.27 0.97A_3_to_10 = fast_mat_exp(mat=A_3, exp=10) #construct A^10
# preview the first 5 rows and 5 columns of A^10 with scientific notation
format(A_3_to_10[1:5,1:5], scientific=TRUE, digits = 3) ## [,1] [,2] [,3] [,4] [,5]
## [1,] "1.35e+12" "1.71e+12" "1.82e+12" "2.22e+12" "1.63e+12"
## [2,] "1.71e+12" "2.16e+12" "2.30e+12" "2.80e+12" "2.06e+12"
## [3,] "1.82e+12" "2.30e+12" "2.46e+12" "2.99e+12" "2.20e+12"
## [4,] "2.22e+12" "2.80e+12" "2.99e+12" "3.64e+12" "2.68e+12"
## [5,] "1.63e+12" "2.06e+12" "2.20e+12" "2.68e+12" "1.97e+12"Again we will test the accuracy of our function to one decimal place.
A_3_to_5_manual = A_3 %*% A_3 %*% A_3 %*% A_3 %*% A_3
A_3_to_10_manual = A_3_to_5_manual %*% A_3_to_5_manual
results_3 = (round(A_3_to_10_manual,1) == round(A_3_to_10,1))
table(results_3)## results_3
## FALSE TRUE
## 13 387Here we see more discrepancies due to floating point errors. However, the errors are likely to be relatively small.