Linear algebra calculations

Install required packages
Matrix algebra
Solving linear equations and calculating the inverse
Norms and distances of matrices and vectors
Trace, determinant, rank, null space, condition index
Matrix decompositions
Orthogonal projections
- Direct implementation of \((X^{t} X)^{-1} X^{t}\)
- Numerically stable implementation using the \(QR\)-decomposition
Detach (automatically) loaded packages (if possible)
Get the article source from GitHub

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

Extract the diagonal of a given matrix.

diag(cov(X))

[1] 43.58333 59.33333

Create a diagonal matrix based on a vector.

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

diag(5, nrow=1)

     [,1]
[1,]    5

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.75   60
[2,]  60.00  178

crossprod(Xdot)

       [,1] [,2]
[1,] 130.75   60
[2,]  60.00  178

(1/(N-1)) * SSP

         [,1]     [,2]
[1,] 43.58333 20.00000
[2,] 20.00000 59.33333

(S <- cov(X))

         [,1]     [,2]
[1,] 43.58333 20.00000
[2,] 20.00000 59.33333

Ds <- diag(1/sqrt(diag(S)))
Ds %*% S %*% Ds

          [,1]      [,2]
[1,] 1.0000000 0.3932968
[2,] 0.3932968 1.0000000

cov2cor(S)

          [,1]      [,2]
[1,] 1.0000000 0.3932968
[2,] 0.3932968 1.0000000

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.5101443 0.6307329
[2,] 0.6631876 0.5872341
[3,] 0.2550722 0.4349882
[4,] 0.4846371 0.2609929

X %*% diag(1/colLens)

          [,1]      [,2]
[1,] 0.5101443 0.6307329
[2,] 0.6631876 0.5872341
[3,] 0.2550722 0.4349882
[4,] 0.4846371 0.2609929

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.69536

sqrt(sum(a1^2))

[1] 9.69536

a2 <- c(6, 9, 10, 8, 7)
A  <- cbind(a1, a2)
sqrt(diag(crossprod(A)))

      a1       a2 
 9.69536 18.16590

sqrt(colSums(A^2))

      a1       a2 
 9.69536 18.16590

norm(A, type="F")

[1] 20.59126

sqrt(crossprod(c(A)))

         [,1]
[1,] 20.59126

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,] 35.62303

dist(B, diag=TRUE, upper=TRUE)

         1        2        3        4
1  0.00000 35.62303 20.80865 30.09983
2 35.62303  0.00000 24.61707 28.58321
3 20.80865 24.61707  0.00000 11.09054
4 30.09983 28.58321 11.09054  0.00000

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]  6.994405e-17  8.366571e-17 -2.466777e-17

Mahalanobis-distance

ideal <- c(1, 2, 3)
y1    <- Y[1, ]
y2    <- Y[2, ]
mat   <- rbind(y1, y2)

mahalanobis(mat, ideal, S)

       y1        y2 
 3.192715 12.257083

Sinv <- solve(S)
t(y1-ideal) %*% Sinv %*% (y1-ideal)

         [,1]
[1,] 3.192715

t(y2-ideal) %*% Sinv %*% (y2-ideal)

         [,1]
[1,] 12.25708

mDist <- mahalanobis(Y, ideal, S)
min(mDist)

[1] 0

(idxMin <- which.min(mDist))

[1] 74

Y[idxMin, ]

[1] 1 2 3

idealM <- t(SsqrtInv %*% (ideal - ctr))
crossprod(Xmt[1, ] - t(idealM))

         [,1]
[1,] 3.192715

crossprod(Xmt[2, ] - t(idealM))

         [,1]
[1,] 12.25708

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))

eigen() decomposition
$values
[1] 9.192582 3.807418

$vectors
           [,1]       [,2]
[1,] -0.9819564  0.1891075
[2,] -0.1891075 -0.9819564

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,] 6.661338e-15    0
[2,] 5.773160e-15    0

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.931935

Xplus <- solve(t(X) %*% X) %*% t(X)
base::norm(X, type="2") * base::norm(Xplus, type="2")

[1] 7.931935

evX <- eigen(t(X) %*% X)$values
sqrt(max(evX) / min(evX[evX >= .Machine$double.eps]))

[1] 7.931935

sqrt(evX / min(evX[evX >= .Machine$double.eps]))

[1] 7.931935 1.502511 1.000000

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.58333  20.00000 -14.00000
[2,]  20.00000  59.33333 -23.33333
[3,] -14.00000 -23.33333  18.66667

eigS <- eigen(S)
G    <- eigS$vectors
D    <- diag(eigS$values)
G %*% D %*% t(G)

          [,1]      [,2]      [,3]
[1,]  43.58333  20.00000 -14.00000
[2,]  20.00000  59.33333 -23.33333
[3,] -14.00000 -23.33333  18.66667

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.372708  1.294488 -1.138518
[2,]  1.294488  7.328159 -1.988899
[3,] -1.138518 -1.988899  3.662611

$Binv
            [,1]        [,2]       [,3]
[1,]  0.16819129 -0.01820349 0.04239706
[2,] -0.01820349  0.16201802 0.08232174
[3,]  0.04239706  0.08232174 0.33091129

$k
[1] 7

$acc
[1] 2.313559e-10

sqrtD <- diag(sqrt(eigS$values))
(A <- G %*% sqrtD %*% t(G))

          [,1]      [,2]      [,3]
[1,]  6.372708  1.294488 -1.138518
[2,]  1.294488  7.328159 -1.988899
[3,] -1.138518 -1.988899  3.662611

A %*% A

          [,1]      [,2]      [,3]
[1,]  43.58333  20.00000 -14.00000
[2,]  20.00000  59.33333 -23.33333
[3,] -14.00000 -23.33333  18.66667

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

N <- eigS$vectors %*% sqrt(diag(eigS$values))
N %*% t(N)

          [,1]      [,2]      [,3]
[1,]  43.58333  20.00000 -14.00000
[2,]  20.00000  59.33333 -23.33333
[3,] -14.00000 -23.33333  18.66667

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:pracma))
try(detach(package:mvtnorm))

Get the article source from GitHub

R markdown - markdown - R code - all posts