Experiments with Random Matrices
Purpose
To explore random matrices
> library(RColorBrewer) > cols <- brewer.pal(4, "Set1") |
Generate a random matrix
> m <- 25
> A <- matrix(data = 0, nrow = m, ncol = m)
> set.seed(12345)
> for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ } |
What do the eigenvalues of a random matrix look like?
> plot(eigen(A)$values, pch = 19, col = cols[1]) |
What happens, say, if you take 100 random matrices and superimpose all their eigenvalues in a single plot? If you do this for m = 8,16, 32, 64,…, what pattern is suggested?
m = 8
> set.seed(12345)
> m <- 8
> eigvals <- matrix(data = 0, nrow = m, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- as.real(eigen(A)$values)
+ }
> cols <- rainbow(100)
> k <- 1
> ran <- range(as.vector(eigvals))
> plot.new()
> for (k in 1:100) {
+ par(new = T)
+ plot(eigvals[, k], pch = 19, col = cols[k], ylim = ran, xlab = "",
+ ylab = "", main = "8 by 8 ")
+ } |
m = 16
> set.seed(12345)
> m <- 16
> eigvals <- matrix(data = 0, nrow = m, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- as.real(eigen(A)$values)
+ }
> cols <- rainbow(100)
> k <- 1
> ran <- range(as.vector(eigvals))
> plot.new()
> for (k in 1:100) {
+ par(new = T)
+ plot(eigvals[, k], pch = 19, col = cols[k], ylim = ran, xlab = "",
+ ylab = "", main = "16 by 16 ")
+ } |
m = 32
> set.seed(12345)
> m <- 32
> eigvals <- matrix(data = 0, nrow = m, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- as.real(eigen(A)$values)
+ }
> cols <- rainbow(100)
> k <- 1
> ran <- range(as.vector(eigvals))
> plot.new()
> for (k in 1:100) {
+ par(new = T)
+ plot(eigvals[, k], pch = 19, col = cols[k], ylim = ran, xlab = "",
+ ylab = "", main = "32 by 32 ")
+ } |
m = 64
> set.seed(12345)
> m <- 64
> eigvals <- matrix(data = 0, nrow = m, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- as.real(eigen(A)$values)
+ }
> cols <- rainbow(100)
> k <- 1
> ran <- range(as.vector(eigvals))
> plot.new()
> for (k in 1:100) {
+ par(new = T)
+ plot(eigvals[, k], pch = 19, col = cols[k], ylim = ran, xlab = "",
+ ylab = "", main = "64 by 64 ")
+ } |
m = 128
> set.seed(12345)
> m <- 128
> eigvals <- matrix(data = 0, nrow = m, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- as.real(eigen(A)$values)
+ }
> cols <- rainbow(100)
> k <- 1
> ran <- range(as.vector(eigvals))
> plot.new()
> for (k in 1:100) {
+ par(new = T)
+ plot(eigvals[, k], pch = 19, col = cols[k], ylim = ran, xlab = "",
+ ylab = "", main = "128 by 128 ")
+ } |
The distribution tends to narrow as increase m
All Eigen values in the same plot
> plot(as.vector(eigvals), pch = 19, col = cols[1], ylim = ran, + xlab = "", ylab = "", main = "128 by 128 ") |
m = 8
> set.seed(12345)
> m <- 8
> eigvals <- matrix(data = 0, nrow = 1, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- max(abs(as.real(eigen(A)$values)))
+ }
> ran <- range(as.vector(eigvals))
> plot(as.vector(eigvals), pch = 19, col = "blue", ylim = ran,
+ xlab = "", ylab = "", main = "8 by 8 ") |
m = 16
> set.seed(12345)
> m <- 16
> eigvals <- matrix(data = 0, nrow = 1, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- max(abs(as.real(eigen(A)$values)))
+ }
> ran <- range(as.vector(eigvals))
> plot(as.vector(eigvals), pch = 19, col = "blue", ylim = ran,
+ xlab = "", ylab = "", main = "16 by 16 ") |
m = 32
> set.seed(12345)
> m <- 32
> eigvals <- matrix(data = 0, nrow = 1, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- max(abs(as.real(eigen(A)$values)))
+ }
> ran <- range(as.vector(eigvals))
> plot(as.vector(eigvals), pch = 19, col = "blue", ylim = ran,
+ xlab = "", ylab = "", main = "32 by 32 ") |
m = 64
> set.seed(12345)
> m <- 64
> eigvals <- matrix(data = 0, nrow = 1, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- max(abs(as.real(eigen(A)$values)))
+ }
> ran <- range(as.vector(eigvals))
> plot(as.vector(eigvals), pch = 19, col = "blue", ylim = ran,
+ xlab = "", ylab = "", main = "64 by 64 ") |
m = 128
> set.seed(12345)
> m <- 128
> eigvals <- matrix(data = 0, nrow = 1, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- max(abs(as.real(eigen(A)$values)))
+ }
> ran <- range(as.vector(eigvals))
> plot(as.vector(eigvals), pch = 19, col = "blue", ylim = ran,
+ xlab = "", ylab = "", main = "128 by 128 ") |
What about condition numbers or more simply, the smallest singular value ?
m = 8
> set.seed(12345)
> m <- 8
> eigvals <- matrix(data = 0, nrow = 1, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- min(abs(as.real(eigen(A)$values)))
+ }
> ran <- range(as.vector(eigvals))
> plot(as.vector(eigvals), pch = 19, col = "blue", ylim = ran,
+ xlab = "", ylab = "", main = "8 by 8 ") |
> hist(eigvals, fill = "blue") |
m = 128
> set.seed(12345)
> m <- 128
> eigvals <- matrix(data = 0, nrow = 1, ncol = 100)
> for (k in 1:100) {
+ A <- matrix(data = 0, nrow = m, ncol = m)
+ for (j in 1:m) {
+ A[, j] <- rnorm(mean = 0, sd = 1/sqrt(m), n = m)
+ }
+ eigvals[, k] <- min(abs(as.real(eigen(A)$values)))
+ }
> ran <- range(as.vector(eigvals))
> plot(as.vector(eigvals), pch = 19, col = "blue", ylim = ran,
+ xlab = "", ylab = "", main = "128 by 128 ") |
> hist(eigvals, fill = "blue") |
