Linear algebra calculations

Install required packages

expm, mvtnorm, pracma

wants <- c("expm", "mvtnorm", "pracma")
has   <- wants %in% rownames(installed.packages())
if(any(!has)) install.packages(wants[!has])

Matrix algebra

Transpose

N  <- 4
Q  <- 2
(X <- matrix(c(20, 26, 10, 19, 29, 27, 20, 12), nrow=N, ncol=Q))
     [,1] [,2]
[1,]   20   29
[2,]   26   27
[3,]   10   20
[4,]   19   12
t(X)
     [,1] [,2] [,3] [,4]
[1,]   20   26   10   19
[2,]   29   27   20   12

Extracting the diagnoal and creating a diagonal matrix

diag(cov(X))
[1] 43.58 59.33
diag(1:3)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    2    0
[3,]    0    0    3
diag(2)
     [,1] [,2]
[1,]    1    0
[2,]    0    1

Multiplication

(Xc <- diag(N) - matrix(rep(1/N, N^2), nrow=N))
      [,1]  [,2]  [,3]  [,4]
[1,]  0.75 -0.25 -0.25 -0.25
[2,] -0.25  0.75 -0.25 -0.25
[3,] -0.25 -0.25  0.75 -0.25
[4,] -0.25 -0.25 -0.25  0.75
(Xdot <- Xc %*% X)
      [,1] [,2]
[1,]  1.25    7
[2,]  7.25    5
[3,] -8.75   -2
[4,]  0.25  -10
(SSP <- t(Xdot) %*% Xdot)
      [,1] [,2]
[1,] 130.8   60
[2,]  60.0  178
crossprod(Xdot)
      [,1] [,2]
[1,] 130.8   60
[2,]  60.0  178
(1/(N-1)) * SSP
      [,1]  [,2]
[1,] 43.58 20.00
[2,] 20.00 59.33
(S <- cov(X))
      [,1]  [,2]
[1,] 43.58 20.00
[2,] 20.00 59.33
Ds <- diag(1/sqrt(diag(S)))
Ds %*% S %*% Ds
       [,1]   [,2]
[1,] 1.0000 0.3933
[2,] 0.3933 1.0000
cov2cor(S)
       [,1]   [,2]
[1,] 1.0000 0.3933
[2,] 0.3933 1.0000
b <- 2
a <- c(-2, 1)
sweep(b*X, 2, a, "+")
     [,1] [,2]
[1,]   38   59
[2,]   50   55
[3,]   18   41
[4,]   36   25
colLens <- sqrt(colSums(X^2))
sweep(X, 2, colLens, "/")
       [,1]   [,2]
[1,] 0.5101 0.6307
[2,] 0.6632 0.5872
[3,] 0.2551 0.4350
[4,] 0.4846 0.2610
X %*% diag(1/colLens)
       [,1]   [,2]
[1,] 0.5101 0.6307
[2,] 0.6632 0.5872
[3,] 0.2551 0.4350
[4,] 0.4846 0.2610

Power

B <- cbind(c(1,1,1), c(0,2,0), c(0,0,2))
B %*% B %*% B
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    7    8    0
[3,]    7    0    8
library(expm)
B %^% 3
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    7    8    0
[3,]    7    0    8

Cross product

a <- c(1, 2, 3)
b <- c(4, 5, 6)
library(pracma)
cross(a, b)
[1] -3  6 -3

Solving linear equations and calculating the inverse

Inverse

Y     <- matrix(c(1, 1, 1, -1), nrow=2)
(Yinv <- solve(Y))
     [,1] [,2]
[1,]  0.5  0.5
[2,]  0.5 -0.5
Y %*% Yinv
     [,1] [,2]
[1,]    1    0
[2,]    0    1

Moore-Penrose generalized inverse

library(MASS)
gInv <- ginv(X)
zapsmall(gInv %*% X)
     [,1] [,2]
[1,]    1    0
[2,]    0    1

Solving linear equations

A  <- matrix(c(9, 1, -5, 0), nrow=2)
b  <- c(5, -3)
(x <- solve(A, b))
[1] -3.0 -6.4
A %*% x
     [,1]
[1,]    5
[2,]   -3

Norms and distances of matrices and vectors

Norm

a1 <- c(3, 4, 1, 8, 2)
sqrt(crossprod(a1))
      [,1]
[1,] 9.695
sqrt(sum(a1^2))
[1] 9.695
a2 <- c(6, 9, 10, 8, 7)
A  <- cbind(a1, a2)
sqrt(diag(crossprod(A)))
    a1     a2 
 9.695 18.166 
sqrt(colSums(A^2))
    a1     a2 
 9.695 18.166 
norm(A, type="F")
[1] 20.59
sqrt(crossprod(c(A)))
      [,1]
[1,] 20.59

Distance

Length of difference vector

set.seed(123)
B <- matrix(sample(-20:20, 12, replace=TRUE), ncol=3)
sqrt(crossprod(B[1, ] - B[2, ]))
      [,1]
[1,] 42.73
dist(B, diag=TRUE, upper=TRUE)
      1     2     3     4
1  0.00 42.73 24.56 25.40
2 42.73  0.00 33.12 35.23
3 24.56 33.12  0.00 32.65
4 25.40 35.23 32.65  0.00

Mahalanobis-transformation

library(mvtnorm)
N     <- 100
mu    <- c(-3, 2, 4)
sigma <- matrix(c(4,2,-3, 2,16,-1, -3,-1,9), byrow=TRUE, ncol=3)
Y     <- round(rmvnorm(N, mean=mu, sigma=sigma))
ctr   <- colMeans(Y)
S     <- cov(Y)
Seig  <- eigen(S)
sqrtD <- sqrt(Seig$values)
SsqrtInv <- Seig$vectors %*% diag(1/sqrtD) %*% t(Seig$vectors)

Xdot  <- sweep(Y, 2, ctr, "-")
Xmt   <- t(SsqrtInv %*% t(Xdot))
zapsmall(cov(Xmt))
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
colMeans(Xmt)
[1] -9.962e-17 -1.671e-17 -5.256e-17

Mahalanobis-distance

ideal <- c(1, 2, 3)
y1    <- Y[1, ]
y2    <- Y[2, ]
mat   <- rbind(y1, y2)
mahalanobis(mat, ideal, S)
   y1    y2 
 5.62 13.17 
Sinv <- solve(S)
t(y1-ideal) %*% Sinv %*% (y1-ideal)
     [,1]
[1,] 5.62
t(y2-ideal) %*% Sinv %*% (y2-ideal)
      [,1]
[1,] 13.17
mDist <- mahalanobis(Y, ideal, S)
min(mDist)
[1] 1.14
(idxMin <- which.min(mDist))
[1] 66
Y[idxMin, ]
[1] -1  1  5
idealM <- t(SsqrtInv %*% (ideal - ctr))
crossprod(Xmt[1, ] - t(idealM))
     [,1]
[1,] 5.62
crossprod(Xmt[2, ] - t(idealM))
      [,1]
[1,] 13.17

Trace, determinant, rank, null space, condition index

Trace

(A <- matrix(c(9, 1, 1, 4), nrow=2))
     [,1] [,2]
[1,]    9    1
[2,]    1    4
sum(diag(A))
[1] 13
sum(diag(t(A) %*% A))
[1] 99
sum(diag(A %*% t(A)))
[1] 99
sum(A^2)
[1] 99

Determinant

det(A)
[1] 35
B <- matrix(c(-3, 4, -1, 7), nrow=2)
all.equal(det(A %*% B), det(A) * det(B))
[1] TRUE
det(diag(1:4))
[1] 24
Ainv <- solve(A)
all.equal(1/det(A), det(Ainv))
[1] TRUE

Rank

qrA <- qr(A)
qrA$rank
[1] 2
(eigA <- eigen(A))
$values
[1] 9.193 3.807

$vectors
        [,1]    [,2]
[1,] -0.9820  0.1891
[2,] -0.1891 -0.9820
zapsmall(eigA$vectors %*% t(eigA$vectors))
     [,1] [,2]
[1,]    1    0
[2,]    0    1
sum(eigA$values)
[1] 13
prod(eigA$values)
[1] 35

Null space (kernel)

library(MASS)
Xnull <- Null(X)
t(X) %*% Xnull
           [,1]      [,2]
[1,] -2.762e-15 1.443e-15
[2,] -2.609e-15 4.441e-16

Condition index

We need to specify base::norm() here because after attaching package expm above, there is another function present with the same name.

X <- matrix(c(20, 26, 10, 19, 29, 27, 20, 12, 17, 23, 27, 25), nrow=4)
kappa(X, exact=TRUE)
[1] 7.932
Xplus <- solve(t(X) %*% X) %*% t(X)
base::norm(X, type="2") * base::norm(Xplus, type="2")
[1] 7.932
evX <- eigen(t(X) %*% X)$values
sqrt(max(evX) / min(evX[evX >= .Machine$double.eps]))
[1] 7.932
sqrt(evX / min(evX[evX >= .Machine$double.eps]))
[1] 7.932 1.503 1.000

Matrix decompositions

Eigenvalues and eigenvectors

X  <- matrix(c(20, 26, 10, 19, 29, 27, 20, 12, 17, 23, 27, 25), nrow=4)
(S <- cov(X))
       [,1]   [,2]   [,3]
[1,]  43.58  20.00 -14.00
[2,]  20.00  59.33 -23.33
[3,] -14.00 -23.33  18.67
eigS <- eigen(S)
G    <- eigS$vectors
D    <- diag(eigS$values)
G %*% D %*% t(G)
       [,1]   [,2]   [,3]
[1,]  43.58  20.00 -14.00
[2,]  20.00  59.33 -23.33
[3,] -14.00 -23.33  18.67

Singular value decomposition

svdX <- svd(X)
all.equal(X, svdX$u %*% diag(svdX$d) %*% t(svdX$v))
[1] TRUE
all.equal(sqrt(eigen(t(X) %*% X)$values), svdX$d)
[1] TRUE

Cholesky decomposition

R <- chol(S)
all.equal(S, t(R) %*% R)
[1] TRUE

\(QR\)-decomposition

We need to specify base::qr.Q() here because after attaching package expm above, there is another function present with the same name.

qrX <- qr(X)
Q   <- base::qr.Q(qrX)
R   <- base::qr.R(qrX)
all.equal(X, Q %*% R)
[1] TRUE

Square-root

library(expm)
sqrtm(S)
$B
       [,1]   [,2]   [,3]
[1,]  6.373  1.294 -1.139
[2,]  1.294  7.328 -1.989
[3,] -1.139 -1.989  3.663

$Binv
        [,1]     [,2]    [,3]
[1,]  0.1682 -0.01820 0.04240
[2,] -0.0182  0.16202 0.08232
[3,]  0.0424  0.08232 0.33091

$k
[1] 7

$acc
[1] 2.313e-10
sqrtD <- diag(sqrt(eigS$values))
(A <- G %*% sqrtD %*% t(G))
       [,1]   [,2]   [,3]
[1,]  6.373  1.294 -1.139
[2,]  1.294  7.328 -1.989
[3,] -1.139 -1.989  3.663
A %*% A
       [,1]   [,2]   [,3]
[1,]  43.58  20.00 -14.00
[2,]  20.00  59.33 -23.33
[3,] -14.00 -23.33  18.67

\(X = N N^{t}\)

N <- eigS$vectors %*% sqrt(diag(eigS$values))
N %*% t(N)
       [,1]   [,2]   [,3]
[1,]  43.58  20.00 -14.00
[2,]  20.00  59.33 -23.33
[3,] -14.00 -23.33  18.67

Orthogonal projections

Direct implementation of \((X^{t} X)^{-1} X^{t}\)

X    <- matrix(c(20, 26, 10, 19, 29, 27, 20, 12, 17, 23, 27, 25), nrow=4)
ones <- rep(1, nrow(X))
P1   <- ones %*% solve(t(ones) %*% ones) %*% t(ones)
P1x  <- P1 %*% X
head(P1x)
      [,1] [,2] [,3]
[1,] 18.75   22   23
[2,] 18.75   22   23
[3,] 18.75   22   23
[4,] 18.75   22   23
a  <- ones / sqrt(crossprod(ones))
P2 <- a %*% t(a)
all.equal(P1, P2)
[1] TRUE
IP1  <- diag(nrow(X)) - P1
IP1x <- IP1 %*% X
all.equal(IP1x, sweep(X, 2, colMeans(X), "-"))
[1] TRUE
A   <- cbind(c(1, 0, 0), c(0, 1, 0))
P3  <- A %*% solve(t(A) %*% A) %*% t(A)
Px3 <- t(P3 %*% t(X))
Px3[1:3, ]
     [,1] [,2] [,3]
[1,]   20   29    0
[2,]   26   27    0
[3,]   10   20    0

Numerically stable implementation using the \(QR\)-decomposition

We need to specify base::qr.Q() here because after attaching package expm above, there is another function present with the same name.

qrX   <- qr(X)
Q     <- base::qr.Q(qrX)
R     <- base::qr.R(qrX)
Xplus <- solve(t(X) %*% X) %*% t(X)
all.equal(Xplus, solve(R) %*% t(Q))
[1] TRUE
all.equal(X %*% Xplus, tcrossprod(Q))
[1] TRUE

Detach (automatically) loaded packages (if possible)

try(detach(package:MASS))
try(detach(package:expm))
try(detach(package:Matrix))
try(detach(package:lattice))
try(detach(package:pracma))
try(detach(package:mvtnorm))

Get the article source from GitHub

R markdown - markdown - R code - all posts