Helper functions for visualization

library(ggplot2)
heatmap = function(A){
  n = nrow(A)
  p = ncol(A)  
  df = data.frame(value = c(A),  i = 1:n, j = rep(1:p, rep(n, p)))
  ggplot(df, aes(j, i, fill = value)) + 
    geom_tile(color = "black")+
    scale_fill_gradient2(low="#619CFF", mid="white", high="#F8766D")+
    scale_y_reverse()+
    theme_void()
}

barplot = function(values){
  n = length(values)
  df = data.frame(value = values,  index = 1:n)
  ggplot(df, aes(index, value, fill = value)) + 
    geom_bar(color = "black", stat = "identity")+
    scale_fill_gradient2(low="#619CFF", mid="white", high="#F8766D")+
    theme_bw()
}

Generating matrices with various structures

Generate a random matrix.

n = 5
p = 3
A = matrix(rnorm(n * p), n, p)
A
##            [,1]       [,2]       [,3]
## [1,] -0.5947032 -1.1770342 -1.3786120
## [2,] -1.6258764  0.2251062 -0.9087142
## [3,]  0.1233574 -0.8423187 -0.2637935
## [4,] -0.5934579  0.1170537  1.4542144
## [5,]  0.5815741  0.7041211 -2.3752939

Find QR decomposition. QR can be used to generate a random column orthogonal matrix.

QR = qr(A)
Q = qr.Q(QR)
Q
##             [,1]       [,2]       [,3]
## [1,] -0.30905523  0.6774535 -0.4059589
## [2,] -0.84493512 -0.2967291 -0.2847709
## [3,]  0.06410637  0.5375551 -0.0675040
## [4,] -0.30840806 -0.1300826  0.4454435
## [5,]  0.30223229 -0.3835644 -0.7423788
t(Q) %*% Q
##               [,1]         [,2]          [,3]
## [1,]  1.000000e+00 8.326673e-17 -8.326673e-17
## [2,]  8.326673e-17 1.000000e+00  5.551115e-17
## [3,] -8.326673e-17 5.551115e-17  1.000000e+00
heatmap(t(Q) %*% Q)

heatmap(Q %*% t(Q))

R = qr.R(QR)
R
##          [,1]       [,2]        [,3]
## [1,] 1.924262  0.2962783  0.01057897
## [2,] 0.000000 -1.6022767 -0.08419703
## [3,] 0.000000  0.0000000  3.24738031
heatmap(R)

Generate a symmetric matrix.

B = matrix(rnorm(n * n), n, n)
B = t(B) + B
heatmap(B)

Compute eigen decomposition.

ED = eigen(B)
lambdas = ED$values
lambdas
## [1]  5.57181240  1.19429049  0.08671586 -1.46284529 -1.65508275
barplot(lambdas)

U = ED$vectors
heatmap(t(U) %*% U)

heatmap(U %*% t(U))

Generate a rank-\(r\) matrix.

r = 1
A = matrix(rnorm(n * r), n, r)
B = matrix(rnorm(p * r), p, r)
C = A %*% t(B)
QR = qr(C)
QR$rank
## [1] 1
heatmap(qr.R(QR))

Generate a PSD matrix.

A = matrix(rnorm(n * p), n, p)
C = A %*% t(A)
heatmap(C)

barplot(eigen(C)$values)

Generate a PD matrix.

A = matrix(rnorm(n * n), n, n)
C = t(A) %*% A
eigen(C)$values
## [1] 13.691791663 10.404435221  5.903383182  4.721992713  0.001418264
barplot(eigen(C)$values)

Find square root of a matrix.

ED = eigen(C)
Csqrt = ED$vectors %*% diag(sqrt(ED$values)) %*% t(ED$vectors)
heatmap(Csqrt %*% Csqrt)

heatmap(C)

Generate a matrix with given eigenvalues.

A = matrix(rnorm(n * n), n, n)
U = qr.Q(qr(A))
lambdas = c(1, 2, 3, 4, 4)
C = U %*% diag(lambdas) %*% t(U)
barplot(eigen(C)$values)

Generate a projection matrix.

A = matrix(rnorm(n * p), n, p)
P = A %*% solve(t(A) %*% A) %*% t(A)
heatmap(P)

barplot(eigen(P)$values)

barplot(diag(P))

Find orthogonal complement to a basis

A = matrix(rnorm(n * p), n, p)
Q = qr.Q(qr(A))
P = Q %*% t(Q)
barplot(eigen(P)$values)

Pperp = diag(n) - P
barplot(eigen(Pperp)$values)

Qperp = eigen(Pperp)$vector[, 1:(n - p)]
t(Qperp) %*% Q
##               [,1]          [,2]          [,3]
## [1,]  1.474515e-17 -1.110223e-16 -5.551115e-17
## [2,] -2.894820e-17 -1.387779e-17  1.249001e-16

Compute singular value decomposition.

A = matrix(rnorm(n * p), n, p)
SVD = svd(A)
U = SVD$u
V = SVD$v
d = SVD$d
heatmap(U)

heatmap(V)

barplot(d)

What if p>n?

A = t(A)
SVD = svd(A)
U = SVD$u
V = SVD$v
d = SVD$d
heatmap(U)

heatmap(V)

barplot(d)

Link between SVD and ED.

A = matrix(rnorm(n * p), n, p)
svd(A)$v
##            [,1]      [,2]       [,3]
## [1,]  0.6810896 0.6906735 -0.2430784
## [2,] -0.1185394 0.4316154  0.8942352
## [3,] -0.7225409 0.5802399 -0.3758408
eigen(t(A) %*% A)$vectors
##            [,1]       [,2]       [,3]
## [1,]  0.6810896 -0.6906735  0.2430784
## [2,] -0.1185394 -0.4316154 -0.8942352
## [3,] -0.7225409 -0.5802399  0.3758408
svd(A)$u
##             [,1]        [,2]        [,3]
## [1,]  0.07988774 -0.44863303  0.02768148
## [2,] -0.40963667  0.46726571 -0.57219423
## [3,] -0.46803574 -0.73316642 -0.16351464
## [4,]  0.50006541 -0.20588474 -0.79120176
## [5,]  0.59723771 -0.02167111  0.13816741
eigen(A %*% t(A))$vectors
##             [,1]        [,2]        [,3]       [,4]        [,5]
## [1,]  0.07988774  0.44863303  0.02768148  0.1123228  0.88258919
## [2,] -0.40963667 -0.46726571 -0.57219423  0.4827467  0.23110581
## [3,] -0.46803574  0.73316642 -0.16351464  0.2921664 -0.36236902
## [4,]  0.50006541  0.20588474 -0.79120176 -0.2707956 -0.09063971
## [5,]  0.59723771  0.02167111  0.13816741  0.7717821 -0.16762925

Centering and computing sample covariance matrix

Matrix view for centering.

X = matrix(rnorm(n * p), n, p)
xbar = colMeans(X)
xbar
## [1]  0.60902145 -0.21215623 -0.08741249
Xc = scale(X, center = T, scale = F)
colMeans(Xc)
## [1] 4.440892e-17 2.775558e-18 0.000000e+00
C = diag(n) - matrix(1/n, n, n)
heatmap(C)

Xc = C %*% X
colMeans(Xc)
## [1] -4.440892e-17  2.081668e-17 -1.665335e-17

Sample covariance matrix.

S = (t(X) %*% C %*% X)/(n - 1)
heatmap(S)

heatmap(cov(X))