Rcode ILA Ch 2

Chapter 2: Matrix Theory and Visualization

Create your working directory named LinAlg and go there

Run the following, replacing “alyss/OneDrive/Documents” with your directory, to set the working directory: setwd(“/Users/alyss/OneDrive/Documents/LinAlg”).

Run getwd() to verify that you are in the directory/folder.

To save a figure

Before plotting the figure, run “setEPS(),” then run “postscript(”figname.eps”, width = w, height = h),” replacing figname, w, and h with the desired parameters. The figure will not show up as an output, and will instead be a file on your computer.

Computational Examples of Matrices

A = matrix(c(1,0,0,4,3, 2), nrow = 3, byrow = TRUE)
B = matrix(c(0,1,-1,2), nrow = 2) #form a matrix by columns
C = A%*%B #matrix multiplication
C

##      [,1] [,2]
## [1,]    0   -1
## [2,]    4    8
## [3,]    2    1

t(C) # transpose matrix of C

##      [,1] [,2] [,3]
## [1,]    0    4    2
## [2,]   -1    8    1

A = matrix(c(1, -1, 1, 2), nrow =2, byrow = TRUE)
A

##      [,1] [,2]
## [1,]    1   -1
## [2,]    1    2

A[,1]

## [1] 1 1

A[1,]

## [1]  1 -1

solve(A) #compute the inverse of A

##            [,1]      [,2]
## [1,]  0.6666667 0.3333333
## [2,] -0.3333333 0.3333333

A%*%solve(A) #verify the inverse of A

##              [,1] [,2]
## [1,] 1.000000e+00    0
## [2,] 1.110223e-16    1

#Solve linear equations
A = matrix(c(30, 40, 1, 1), 
           nrow =2, byrow = TRUE)
b = c(1000, 30)
solve(A,b)

## [1] 20 10

solve(A)%*%b #Another way to solve the equations

##      [,1]
## [1,]   20
## [2,]   10

det(A) #compute the determinant

## [1] -10

#Dot product of two vectors
u = c(2,0,1)
v = c(2,-4,3)
t(u)%*%v

##      [,1]
## [1,]    7

library(pracma)

## Warning: package 'pracma' was built under R version 4.4.3

dot(u,v) #for nD

## [1] 7

cross(u,v) #only for 3D

## [1]  4 -4 -8

library(Matrix)

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:pracma':
## 
##     expm, lu, tril, triu

rankMatrix(A) #Find the rank of a matrix

## [1] 2
## attr(,"method")
## [1] "tolNorm2"
## attr(,"useGrad")
## [1] FALSE
## attr(,"tol")
## [1] 4.440892e-16

#Orthogonal matrices
p = sqrt(2)/2
Q = matrix(c(p,-p,p,p), nrow=2) 
Q #is an orthogonal matrix 

##            [,1]      [,2]
## [1,]  0.7071068 0.7071068
## [2,] -0.7071068 0.7071068

Q%*%t(Q) #verify O as an orthogonal matrix

##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

det(Q) #The determinant of an orthogonal matrix is 1 or -1

## [1] 1

#Generate a 2-by-4 zero matrix 
zeroM = matrix(0, nrow = 2, ncol = 4)
zeroM

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]    0    0    0    0

matrix(rep(0, 8), nrow =2)

##      [,1] [,2] [,3] [,4]
## [1,]    0    0    0    0
## [2,]    0    0    0    0

#Generate a 3-by-3 identity matrix
diag(1, 3)

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

diag(3)

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

#generate a 3-by-3 diagonal matrix
diag(c(2.4, 3.1, -0.9))

##      [,1] [,2] [,3]
## [1,]  2.4  0.0  0.0
## [2,]  0.0  3.1  0.0
## [3,]  0.0  0.0 -0.9

Eigenvalues and Eigenvectors

A = matrix(c(1, 2, 2, 1), nrow=2)
eigen(A)

## eigen() decomposition
## $values
## [1]  3 -1
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068

#R plot eigenvector v vs a non-eigenvector u
setwd("/Users/alyss/OneDrive/Documents/LinAlg")

par(mar=c(4.5,4.5,2.0,0.5))
plot(9,9,
     main = 'An eigenvector vs a non-eigenvector',
     cex.axis = 1.4, cex.lab = 1.4,
     xlim = c(0,3), ylim=c(0,3),
     xlab = bquote(x[1]), ylab = bquote(x[2]))
arrows(0,0, 1,0, length = 0.25, 
       angle = 8, lwd = 5, col = 'blue')
arrows(0,0, 1,2, length = 0.3, 
       angle = 8, lwd = 2, col = 'blue',  lty = 3)
arrows(0,0, 1,1, length = 0.25, 
       angle = 8, lwd = 5, col='red') 
arrows(0,0, 3,3, length = 0.3, 
       angle = 8, lwd = 2, col='red', lty = 3)
text(1.4,0.1, 'Non-eigenvector u', cex =1.4, col = 'blue')
text(1.0,2.1, 'Au', cex =1.4, col = 'blue')
text(1.5,0.9, 'Eigenvector v', cex =1.4, col = 'red')
text(2.8, 2.95, 'Av', cex =1.4, col = 'red')

#dev.off()

# Verify diagonalization and decomposition: R code
C = matrix(c(2,1,1,2), nrow = 2)
eigen(C)

## eigen() decomposition
## $values
## [1] 3 1
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068

Q = eigen(C)$vectors
D = t(Q)%*%C%*%Q #Matrix diagonalization
D

##      [,1] [,2]
## [1,]    3    0
## [2,]    0    1

Q%*%D%*%t(Q) #Matrix decomposition

##      [,1] [,2]
## [1,]    2    1
## [2,]    1    2

D[1,1]*Q[,1]%*%t(Q[,1]) + D[2,2]*Q[,2]%*%t(Q[,2])

##      [,1] [,2]
## [1,]    2    1
## [2,]    1    2

Some Types of Matrix and Vector Products

#Hadamard product of two matrices
#install.packages('matrixcalc')
library(matrixcalc)
A = matrix( c( 1, 2, 3, 4 ), nrow=2, byrow=TRUE )
B = matrix( c( 2, 4, 6, 8 ), nrow=2, byrow=TRUE )
hadamard.prod(A, B)

##      [,1] [,2]
## [1,]    2    8
## [2,]   18   32

#Jordan product of A and B
A = matrix( c( 1, 2, 3, 4 ), nrow=2, byrow=TRUE )
B = matrix( c( 2, 1, 2, 1 ), nrow=2, byrow=TRUE )
(A%*%B + B%*%A)/2

##      [,1] [,2]
## [1,]  5.5  5.5
## [2,]  9.5  7.5

#R commutator of A and B
A = matrix( c( 1, 2, 3, 4 ), nrow=2, byrow=TRUE )
B = matrix( c( 2, 1, 2, 1 ), nrow=2, byrow=TRUE )
A%*%B - B%*%A

##      [,1] [,2]
## [1,]    1   -5
## [2,]    9   -1

#install.packages('psych')
library(psych)

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:pracma':
## 
##     logit, polar

tr(A%*%B - B%*%A) #tr for trace

## [1] 0

#Cross product
library(pracma)
x = 1:3
y = 4:6
cross(x, y) 

## [1] -3  6 -3

dot(x, y)

## [1] 32

#Outer product of two vectors
a = 1:2
b = 1:4
a%o%b #outer product a_2-by-1 times t(b)_1-by-4

##      [,1] [,2] [,3] [,4]
## [1,]    1    2    3    4
## [2,]    2    4    6    8

#Outer product of A_mn and B_nq
A = matrix(1:4, ncol = 2)
B = matrix(1:6, ncol = 3)
A%o%B

## , , 1, 1
## 
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4
## 
## , , 2, 1
## 
##      [,1] [,2]
## [1,]    2    6
## [2,]    4    8
## 
## , , 1, 2
## 
##      [,1] [,2]
## [1,]    3    9
## [2,]    6   12
## 
## , , 2, 2
## 
##      [,1] [,2]
## [1,]    4   12
## [2,]    8   16
## 
## , , 1, 3
## 
##      [,1] [,2]
## [1,]    5   15
## [2,]   10   20
## 
## , , 2, 3
## 
##      [,1] [,2]
## [1,]    6   18
## [2,]   12   24

dim(A%o%B)

## [1] 2 2 2 3

#Outer product of A_mn and B_pq
A = matrix(1:4, ncol = 2)
B = matrix(1:9, ncol = 3)
A%o%B

## , , 1, 1
## 
##      [,1] [,2]
## [1,]    1    3
## [2,]    2    4
## 
## , , 2, 1
## 
##      [,1] [,2]
## [1,]    2    6
## [2,]    4    8
## 
## , , 3, 1
## 
##      [,1] [,2]
## [1,]    3    9
## [2,]    6   12
## 
## , , 1, 2
## 
##      [,1] [,2]
## [1,]    4   12
## [2,]    8   16
## 
## , , 2, 2
## 
##      [,1] [,2]
## [1,]    5   15
## [2,]   10   20
## 
## , , 3, 2
## 
##      [,1] [,2]
## [1,]    6   18
## [2,]   12   24
## 
## , , 1, 3
## 
##      [,1] [,2]
## [1,]    7   21
## [2,]   14   28
## 
## , , 2, 3
## 
##      [,1] [,2]
## [1,]    8   24
## [2,]   16   32
## 
## , , 3, 3
## 
##      [,1] [,2]
## [1,]    9   27
## [2,]   18   36

dim(A%o%B)

## [1] 2 2 3 3

#Kronecker product
library(fastmatrix)

## 
## Attaching package: 'fastmatrix'

## The following object is masked from 'package:psych':
## 
##     minkowski

## The following objects are masked from 'package:matrixcalc':
## 
##     vec, vech

## The following object is masked from 'package:Matrix':
## 
##     lu

## The following objects are masked from 'package:pracma':
## 
##     geomean, hadamard, lu

A <- diag(1:2)
B <- matrix(1:4, ncol = 2)
kronecker.prod(A, B)

##      [,1] [,2] [,3] [,4]
## [1,]    1    3    0    0
## [2,]    2    4    0    0
## [3,]    0    0    2    6
## [4,]    0    0    4    8

# an example with vectors
ones <- rep(1, 2)
y <- 1:4
kronecker.prod(ones, t(y)) # 2-by-4 matrix

##      [,1] [,2] [,3] [,4]
## [1,]    1    2    3    4
## [2,]    1    2    3    4

#Cholesky decomposition
dat = matrix(1:4, ncol = 2)
A = dat%*%t(dat)
chol(A)

##          [,1]      [,2]
## [1,] 3.162278 4.4271887
## [2,] 0.000000 0.6324555

#Direct sum of two matrices
A = matrix(1:4, ncol = 2)
B = matrix(1:6, ncol = 3)
A1 = rbind(A, matrix(rep(0, 4), ncol = 2))
B1 = rbind(matrix(rep(0, 6), ncol = 3), B)
C = cbind(A1, B1) #= direct sum of A and B
C

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    3    0    0    0
## [2,]    2    4    0    0    0
## [3,]    0    0    1    3    5
## [4,]    0    0    2    4    6

#Express the direct sum by Kronecker products
#kronecker.prod(diag(1,0), A) + kronecker.prod(diag(0,1), B)

SVD Decompoisitions and Reconstructions

#SVD example for a 2-by-3 matrix
A=matrix(c(-1,1,0,2,-2,3),nrow=2)
A #Show the 2-by-3 matrix

##      [,1] [,2] [,3]
## [1,]   -1    0   -2
## [2,]    1    2    3

svdA=svd(A) #Compute the SVD of A and put the results in svdA
svdA #Show SVD results: d, U, and V

## $d
## [1] 4.221571 1.085514
## 
## $u
##            [,1]      [,2]
## [1,] -0.4791881 0.8777123
## [2,]  0.8777123 0.4791881
## 
## $v
##           [,1]       [,2]
## [1,] 0.3214207 -0.3671293
## [2,] 0.4158226  0.8828774
## [3,] 0.8507528 -0.2928200

round(svdA$d, digits=2) #Show only the singular values

## [1] 4.22 1.09

round(svdA$u, digits=2) #Show only matrix U

##       [,1] [,2]
## [1,] -0.48 0.88
## [2,]  0.88 0.48

round(svdA$v, digits=2)#Show only matrix V

##      [,1]  [,2]
## [1,] 0.32 -0.37
## [2,] 0.42  0.88
## [3,] 0.85 -0.29

sqrt(eigen(A%*%t(A))$values)

## [1] 4.221571 1.085514

#Data reconstruction by singular vectors: R code
round(svdA$d[1]*svdA$u[,1]%*%t(svdA$v[,1]), digits=1)

##      [,1] [,2] [,3]
## [1,] -0.7 -0.8 -1.7
## [2,]  1.2  1.5  3.2

round(svdA$d[1]*svdA$u[,1]%*%t(svdA$v[,1]) + 
        svdA$d[2]*svdA$u[,2]%*%t(svdA$v[,2]), 
      digits =2)

##      [,1] [,2] [,3]
## [1,]   -1    0   -2
## [2,]    1    2    3

#R plot schematic diagram of SVD
setwd("/Users/alyss/OneDrive/Documents/LinAlg")

par(mar=c(0,0,0,0))
plot(200, axes = FALSE,
     xlab = "", ylab = "",
     xlim = c(-3,28), ylim = c(-3,16))
text(13,15.5, cex=2.2,
     bquote("SVD:" ~ A==UDV^t~ "when n > m or n < m"))
#Space-time data matrix A when n>m
segments(x0 = c(0,0,3,3),
         y0 = c(6,12,12,6) +1,
         x1 = c(0,3,3,0),
         y1 = c(12,12,6,6) +1, 
         col = c('blue','red','blue','red'),lwd =3)
segments(x0 = c(0.5,1.0),
         y0 = c(6,6)+1,
         x1 = c(0.5,1.0),
         y1 = c(12,12)+1,
         lwd =1.3, lty = 3)
text(-.8, 9+1, 'n', srt=90, col ='blue', cex = 1.4)
text(1.5, 12.8+1, 'm', col = 'red',  cex = 1.4)
text(2.0, 9+1, '...',  cex = 1.4)
text(2, 5+1, bquote(A[n%*%m]),  cex = 2.5)
text(5, 9+1, '=',  cex = 3)
#Spatial matrix U
segments(x0 = c(7,7,10,10),
         y0 = c(6,12,12,6)+1,
         x1 = c(7,10,10,7),
         y1 = c(12,12,6,6)+1, 
         col = c('blue','blue','blue','blue'), lwd =3)
segments(x0 = c(7.5,8),
         y0 = c(6,6)+1,
         x1 = c(7.5,8),
         y1 = c(12,12)+1,
         lwd =1.3, lty = 3, col = 'blue')
text(6.2, 9+1, 'n', srt=90, col ='blue', cex = 1.4)
text(8.5, 12.8+1, 'm', col = 'red',  cex = 1.4)
text(9, 9+1, '...',  cex = 1.4, col='blue')
text(8.7, 5.0+1, bquote(U[n%*%m]),  cex = 2.5, col= 'blue')
#Singular value diagonal matrix D
segments(x0 = c(12,12,15,15),
         y0 = c(9,12,12,9)+1,
         x1 = c(12,15,15,12),
         y1 = c(12,12,9,9)+1, 
         col = c('brown','brown','brown','brown'), lwd =3)
segments(x0 = 12, y0 = 12+1, x1 = 15, y1 = 9+1, lty=3,
         col = c('brown'), lwd =1.3)#diagonal line
text(11.2, 10.5+1, 'm', srt=90, col ='red', cex = 1.4)
text(13.5, 12.8+1, 'm', col = 'red',  cex = 1.4)
text(14.1, 11.3+1, '0', col = 'brown',  cex = 1.4)
text(12.9, 10.0+1, '0', col = 'brown',  cex = 1.4)
text(13.9, 8.0+1, bquote(D[m%*%m]),  cex = 2.5, col='brown')
#Temporal matrix V
segments(x0 = c(17,17,20,20),
         y0 = c(9,12,12,9)+1,
         x1 = c(17,20,20,17),
         y1 = c(12,12,9,9)+1, 
         col = c('red','red','red','red'), lwd =3)
segments(x0 = c(17,17),
         y0 = c(11.5,10.8)+1,
         x1 = c(20,20),
         y1 = c(11.5,10.8)+1, 
         col = c('red','red'), lty=3, lwd =1.3)
text(16.2, 10.5+1, 'm', srt=90, col ='red', cex = 1.4)
text(18.5, 12.5+1, 'm', col = 'red',  cex = 1.4)
text(19.5, 8+1, bquote((V^t)[m%*%m]),  cex = 2.5, col='red')
text(18.5, 10+1, '...',  col='red', srt=90, cex =1.4)
# Space-time data matrix B when n < m
segments(x0 = c(0,0,6,6),
         y0 = c(0,3,3,0),
         x1 = c(0,6,6,0),
         y1 = c(3,3,0,0), 
         col = c('blue','red','blue','red'), lwd =3)
segments(x0 = c(1,2,5),
         y0 = c(0,0,0),
         x1 = c(1,2,5),
         y1 = c(3,3,3),
         lwd =1.3, lty = 3)
text(-0.8, 1.5, 'n', srt=90, col ='blue', cex = 1.4)
text(3, 3.8, 'm', col = 'red',  cex = 1.4)
text(3.5, 1.5, '...',  cex = 1.4)
text(3, -1.5, bquote(A[n%*%m]),  cex = 2.5)
text(8, 1.5, '=',  cex = 3)
#Spatial matrix U
segments(x0 = c(11,11,14,14),
         y0 = c(0,3,3,0),
         x1 = c(11,14,14,11),
         y1 = c(3,3,0,0), 
         col = c('blue','blue','blue','blue'), lwd =3)
segments(x0 = c(11.5,12.2),
         y0 = c(0,0),
         x1 = c(11.5,12.2),
         y1 = c(3,3),
         lwd =1.3, lty = 3, col = 'blue')
text(10.2, 1.5, 'n', srt=90, col ='blue', cex = 1.4)
text(12.5, 3.8, 'n', col = 'blue',  cex = 1.4)
text(13.2, 1.5, '...',  cex = 1.4, col='blue')
text(12.5, -1.5, bquote(U[n%*%n]),  cex = 2.5, col= 'blue')
#Singular value diagonal matrix D
segments(x0 = c(16,16,19,19),
         y0 = c(0,3,3,0),
         x1 = c(16,19,19,16),
         y1 = c(3,3,0,0), 
         col = c('brown','brown','brown','brown'), lwd =3)
segments(x0 = 16, y0 = 3, x1 = 19, y1 = 0, lty=3,
         col = c('brown'), lwd =1.3)#diagonal line
text(15.2, 1.5, 'n', srt=90, col ='blue', cex = 1.4)
text(17.5, 3.8, 'n', col = 'blue',  cex = 1.4)
text(18.1, 2.3, '0', col = 'brown',  cex = 1.4)
text(16.9, 1.0, '0', col = 'brown',  cex = 1.4)
text(17.5, -1.5, bquote(D[n%*%n]),  cex = 2.5, col='brown')
#Temporal matrix V
segments(x0 = c(21,21,27,27),
         y0 = c(0,3,3,0),
         x1 = c(21,27,27,21),
         y1 = c(3,3,0,0), 
         col = c('red','red','red','red'),
         lwd =3)
segments(x0 = c(21,21),
         y0 = c(2.5,1.8),
         x1 = c(27,27),
         y1 = c(2.5,1.8), 
         col = c('red','red'), lty=3, lwd =1.3)
text(20.2, 1.5, 'n', srt=90, col ='blue', cex = 1.4)
text(24, 3.8, 'm', col = 'red',  cex = 1.4)
text(24, -1.5, bquote((V^t)[n%*%m]),  cex = 2.5, col='red')
text(24, 1, '...',  col='red', srt=90, cex =1.4)

#dev.off()

SVD Example: Darwin/Tahiti Data

#R SVD analysis for the weighted SOI from SLP data
#for seven years: 2009-2015
setwd("/Users/alyss/OneDrive/Documents/LinAlg")
Pda<-read.table("data/PSTANDdarwin.txt", header=F)
dim(Pda) #Monthly Darwin data from 1951-2015

## [1] 65 13

Pda

##      V1   V2   V3   V4   V5   V6   V7   V8   V9  V10  V11  V12  V13
## 1  1951 -1.3 -1.6 -1.2 -0.4  0.4 -1.6 -0.1 -0.2  0.2  0.6  0.5  1.4
## 2  1952  0.0  0.2  0.0  0.4 -1.1 -0.4 -0.4 -0.7 -0.1 -0.8 -1.1  0.9
## 3  1953  0.2  0.4 -0.3 -0.9  1.5 -0.4  0.4 -0.3  0.9  0.0  0.5  0.6
## 4  1954 -1.6  0.4 -0.6 -2.0 -1.3 -0.6 -0.9 -1.4 -0.3 -0.1  0.2 -0.8
## 5  1955  1.5 -2.8 -1.5 -0.9 -1.9 -1.8 -1.6 -1.8 -1.2 -1.3 -1.2  0.2
## 6  1956 -1.7 -1.9 -1.3 -1.4 -2.2 -1.5 -2.0 -1.7 -0.8 -1.0 -0.2 -1.6
## 7  1957 -0.2 -0.4 -0.5 -0.4  0.7 -0.4 -0.5 -0.1  0.8  0.7  0.8  0.3
## 8  1958  1.7  0.2 -0.2 -0.9 -0.5 -1.2 -0.9 -1.6 -0.7 -1.1  0.2  0.8
## 9  1959 -0.5  1.7 -1.3 -0.9 -1.0 -0.1  0.1 -0.3 -0.3 -0.2 -1.1 -1.2
## 10 1960 -0.2 -0.1 -1.1 -1.3 -1.4 -0.4 -0.8 -0.8 -1.6 -0.1 -0.2 -0.2
## 11 1961  0.5 -0.4  1.1 -0.6 -0.4 -0.1 -0.3 -0.2 -0.1 -0.3  0.3 -0.6
## 12 1962 -1.4  0.5  0.1 -0.1 -1.3 -1.3 -1.0 -0.1 -1.0 -1.0 -0.7 -0.1
## 13 1963 -0.7  1.1 -1.5 -1.3 -1.1  0.3  0.1 -0.9  0.4  1.5  1.3  1.5
## 14 1964  0.6  0.2 -1.2 -0.7 -0.9 -1.0 -1.3 -1.2 -1.7 -1.1 -1.0  0.2
## 15 1965 -0.3 -1.0 -1.3  1.2 -0.3  1.0  1.9  0.3  1.0  1.5  1.4 -0.7
## 16 1966  0.7 -0.1  0.2  0.0  0.7  0.3 -0.5 -0.7 -0.1  0.5  0.0  0.2
## 17 1967 -1.5 -1.1 -1.0  0.7 -0.4 -0.4 -0.6 -0.5  0.4  0.0  0.0  0.4
## 18 1968 -2.3 -1.6 -0.1 -0.2 -1.8 -0.7 -0.7 -0.6  0.0  0.2  0.8 -0.1
## 19 1969 -0.8 -0.5  0.0  0.9 -0.4 -0.2  0.1  0.0 -0.4  0.7  0.2 -0.4
## 20 1970  1.4  0.3 -1.0  0.1 -0.8 -1.0  0.3 -1.3 -1.9 -1.9 -1.3 -2.0
## 21 1971 -0.8 -2.1 -2.1 -1.7 -1.0  0.1 -0.8 -1.5 -1.3 -1.1 -0.5 -1.4
## 22 1972 -0.8 -1.6 -1.5  0.6  2.0 -0.1  0.3  0.4  1.3  1.7  0.8  2.0
## 23 1973  0.6  1.3 -1.2 -0.8 -0.6 -1.8 -0.9 -1.7 -0.8 -0.7 -2.9 -0.7
## 24 1974 -2.4 -2.2 -2.6 -1.0 -0.4 -0.5 -0.8 -1.2 -1.1 -1.1 -1.1 -1.2
## 25 1975 -0.3  0.1 -1.6 -1.9 -1.0 -1.4 -1.6 -1.9 -1.7 -1.7 -0.8 -1.7
## 26 1976 -1.4 -1.3 -2.5  0.0  0.5 -0.3  0.5  0.7  0.6 -0.8 -1.0 -0.8
## 27 1977 -0.7 -1.1 -0.2 -0.5 -0.2  0.5  0.2  0.2  0.6  1.4  1.0  0.6
## 28 1978  1.4  1.9  1.1  0.2 -2.0 -1.7 -1.9 -1.3 -0.6 -0.1  0.2  0.0
## 29 1979 -0.6  0.1 -0.4 -0.2 -0.3 -0.3 -2.1 -0.1 -0.5  0.2  0.2  1.0
## 30 1980 -0.7 -0.5 -1.0 -0.2 -0.2 -0.4 -0.2  0.0  0.4 -0.1  0.2 -0.1
## 31 1981 -1.4 -0.6  1.4  0.0 -0.9 -1.3 -1.5 -0.9 -0.7  0.3 -0.5  0.1
## 32 1982 -0.3  0.2 -0.8  0.2  0.4  0.9  1.0  1.2  1.6  1.6  2.0  1.7
## 33 1983  2.8  3.1  1.6  0.2  0.2  0.3  0.9  0.7 -0.1 -0.2  0.7  1.0
## 34 1984  0.4 -1.4 -0.2  0.2 -0.4  0.2 -0.2 -1.1 -0.4  0.5  0.4 -2.1
## 35 1985  0.7 -1.3 -2.0 -2.2 -1.9 -0.3  0.2 -0.9 -0.7 -0.1 -0.5  0.2
## 36 1986 -1.8  0.1  0.1 -0.5  0.1 -1.8 -0.6 -1.3  0.0  0.1  0.4  2.1
## 37 1987  0.4  0.5  2.1  1.5  1.0  0.2  1.1  0.8  1.1  0.4  0.0  0.4
## 38 1988  1.0 -0.1 -1.4  0.2 -0.6  0.9 -0.6 -0.8 -1.5 -1.5 -1.4 -0.8
## 39 1989 -1.2  0.0 -0.9 -2.1 -1.3 -0.1 -1.2 -0.1 -0.5 -0.8 -0.2  1.4
## 40 1990  0.2  2.9  0.0 -0.2 -1.1 -0.5 -0.6  0.1  1.0  0.1  0.2  0.3
## 41 1991 -0.3 -0.3  1.1  0.0  1.0  0.3  0.0  0.6  1.0  0.3  1.1  1.0
## 42 1992  3.1  0.8  2.0  1.0 -0.4 -0.4  1.3 -0.5 -0.9  0.4  0.5 -0.1
## 43 1993  0.1 -0.4  1.6  1.4  0.8  0.9  0.1  1.6  1.1  1.1  0.1 -0.4
## 44 1994  0.0 -0.3  1.0  0.8  0.5  0.0  0.6  1.4  1.1  1.0  0.8  1.7
## 45 1995  0.3  0.9 -1.0 -0.3  0.1  0.1 -0.7 -0.2 -0.2 -0.5 -0.1  0.1
## 46 1996 -0.9 -0.6 -0.8 -0.8  0.6 -0.5 -0.8 -0.9 -1.8 -1.0 -0.5 -1.6
## 47 1997  0.6 -1.8  0.2  1.8  0.4  1.6  0.8  1.6  2.2  2.4  2.3  1.1
## 48 1998  1.0  3.1  1.5  1.1  0.0 -0.8 -0.8 -0.6 -0.6 -1.1 -1.5 -1.4
## 49 1999 -1.6 -0.6 -1.8 -1.1  0.6  0.2 -0.2  0.6  0.1 -0.4 -0.8 -0.7
## 50 2000  0.0 -0.2 -0.8 -1.2  0.2  0.5 -0.3 -1.1 -0.2 -1.5 -2.2 -2.3
## 51 2001  0.5 -3.0 -0.4 -0.4  0.6 -0.7  0.0 -0.1 -0.4 -0.8 -0.5  0.3
## 52 2002  0.5 -0.2  0.7  0.1  0.6  0.2  1.1  0.3  0.7 -0.2  0.8  1.5
## 53 2003  0.2 -0.2  0.1  0.6  0.0  0.1  0.0 -0.3 -0.1 -0.4  0.4 -1.1
## 54 2004  0.8 -0.6 -1.4  0.6 -0.8  1.4  0.4  0.1  0.8  0.3  0.3 -0.1
## 55 2005 -0.3  2.0  0.1  1.1  0.7 -0.8 -0.2  0.1 -0.3 -0.8 -0.5  0.3
## 56 2006 -1.8  0.5 -1.9 -2.0  0.6  0.5  0.9  1.1  1.0  1.9  0.8  1.2
## 57 2007  0.0  0.2 -0.6  0.6  0.3 -1.3  0.7  0.0 -0.7 -0.7 -1.2 -1.3
## 58 2008 -1.6 -1.9  0.1  0.0  0.6  0.2  0.0  0.1 -0.5  0.2 -0.7 -0.8
## 59 2009 -0.6 -1.5  0.1 -0.5 -0.8 -0.2 -0.3 -0.4 -0.8  0.6  0.0  0.5
## 60 2010 -0.7  0.9  0.5 -0.3 -1.2  0.2 -0.9 -1.1 -1.3 -1.3 -0.4 -2.3
## 61 2011 -1.4 -1.7 -1.8 -0.9  0.5  0.5 -0.2  0.2  0.3 -0.4 -0.3 -2.2
## 62 2012 -0.6  0.3 -1.8  0.7 -0.2  0.3 -0.6  0.6  0.0  0.0  0.4  0.3
## 63 2013 -0.8  0.3 -1.0 -0.2 -0.6 -1.6 -0.8 -0.3 -0.6 -0.1 -1.1  0.3
## 64 2014 -1.5 -0.8  0.7 -0.9  0.2 -0.4  0.4  1.0  0.8  0.8  0.9  0.1
## 65 2015 -0.2  0.9  0.5  0.6  1.5  0.3  1.4  1.3  1.5  2.8  0.5 -0.4

pdaDec<-Pda[,13] #Darwin Dec standardized SLP anomalies data
pdaDec

##  [1]  1.4  0.9  0.6 -0.8  0.2 -1.6  0.3  0.8 -1.2 -0.2 -0.6 -0.1  1.5  0.2 -0.7
## [16]  0.2  0.4 -0.1 -0.4 -2.0 -1.4  2.0 -0.7 -1.2 -1.7 -0.8  0.6  0.0  1.0 -0.1
## [31]  0.1  1.7  1.0 -2.1  0.2  2.1  0.4 -0.8  1.4  0.3  1.0 -0.1 -0.4  1.7  0.1
## [46] -1.6  1.1 -1.4 -0.7 -2.3  0.3  1.5 -1.1 -0.1  0.3  1.2 -1.3 -0.8  0.5 -2.3
## [61] -2.2  0.3  0.3  0.1 -0.4

Pta<-read.table("data/PSTANDtahiti.txt", header=F)
ptaDec=Pta[,13] #Tahiti Dec standardized SLP anomalies
ptada1 = cbind(pdaDec, ptaDec) #space-time data matrix
ptada1

##       pdaDec ptaDec
##  [1,]    1.4    0.1
##  [2,]    0.9   -1.1
##  [3,]    0.6   -0.1
##  [4,]   -0.8    1.5
##  [5,]    0.2    1.8
##  [6,]   -1.6    0.1
##  [7,]    0.3   -0.2
##  [8,]    0.8   -0.2
##  [9,]   -1.2    0.3
## [10,]   -0.2    1.1
## [11,]   -0.6    1.8
## [12,]   -0.1    0.2
## [13,]    1.5   -0.6
## [14,]    0.2   -0.3
## [15,]   -0.7   -0.4
## [16,]    0.2   -0.4
## [17,]    0.4   -0.6
## [18,]   -0.1    0.2
## [19,]   -0.4    0.2
## [20,]   -2.0    1.2
## [21,]   -1.4   -1.0
## [22,]    2.0   -0.1
## [23,]   -0.7    2.3
## [24,]   -1.2   -0.8
## [25,]   -1.7    1.7
## [26,]   -0.8   -1.4
## [27,]    0.6   -1.1
## [28,]    0.0   -0.1
## [29,]    1.0   -0.2
## [30,]   -0.1   -0.2
## [31,]    0.1    0.9
## [32,]    1.7   -2.0
## [33,]    1.0    1.0
## [34,]   -2.1   -2.2
## [35,]    0.2    0.6
## [36,]    2.1   -0.3
## [37,]    0.4   -0.4
## [38,]   -0.8    1.2
## [39,]    1.4    0.5
## [40,]    0.3   -0.1
## [41,]    1.0   -1.9
## [42,]   -0.1   -1.0
## [43,]   -0.4    0.0
## [44,]    1.7   -0.3
## [45,]    0.1   -0.7
## [46,]   -1.6   -0.1
## [47,]    1.1   -0.5
## [48,]   -1.4    0.9
## [49,]   -0.7    1.6
## [50,]   -2.3   -0.9
## [51,]    0.3   -1.1
## [52,]    1.5   -0.3
## [53,]   -1.1    0.8
## [54,]   -0.1   -1.4
## [55,]    0.3    0.2
## [56,]    1.2    0.7
## [57,]   -1.3    1.5
## [58,]   -0.8    1.6
## [59,]    0.5   -0.7
## [60,]   -2.3    2.5
## [61,]   -2.2    1.9
## [62,]    0.3   -0.7
## [63,]    0.3    0.4
## [64,]    0.1   -0.8
## [65,]   -0.4   -1.3

#Space-time data format
ptada = t(ptada1[59:65,]) #2009-2015 data
colnames(ptada)<-2009:2015
rownames(ptada)<-c("Darwin", "Tahiti")
ptada #6 year of data for two stations

##        2009 2010 2011 2012 2013 2014 2015
## Darwin  0.5 -2.3 -2.2  0.3  0.3  0.1 -0.4
## Tahiti -0.7  2.5  1.9 -0.7  0.4 -0.8 -1.3

dim(ptada)

## [1] 2 7

svdptd = svd(ptada) #SVD for the 2-by-6 matrix
U=round(svdptd$u, digits=2)
U

##       [,1] [,2]
## [1,] -0.66 0.75
## [2,]  0.75 0.66

D=round(diag(svdptd$d), digits=2)
D

##      [,1] [,2]
## [1,]  4.7 0.00
## [2,]  0.0 1.42

V =round(svdptd$v, digits=2)
t(V)

##       [,1]  [,2]  [,3]  [,4] [,5]  [,6]  [,7]
## [1,] -0.18  0.72  0.61 -0.15 0.02 -0.14 -0.15
## [2,] -0.06 -0.06 -0.28 -0.17 0.34 -0.32 -0.82

#For the entire period of 1951-2015
setwd("/Users/alyss/OneDrive/Documents/LinAlg")
Pda<-read.table("data/PSTANDdarwin.txt", header=F)
dim(Pda) #Monthly Darwin data from 1951-2015

## [1] 65 13

pdaDec<-Pda[,13] #Darwin Dec standardized SLP anomalies data
Pta<-read.table("data/PSTANDtahiti.txt", header=F)
ptaDec=Pta[,13] #Tahiti Dec standardized SLP anomalies
ptada1 = cbind(pdaDec, ptaDec) #space-time data matrix

#Space-time data format
ptada = t(ptada1) #2009-2015 data
colnames(ptada)<-1951:2015
rownames(ptada)<-c("Darwin", "Tahiti")
ptada[,1:6] #6 year of data for two stations

##        1951 1952 1953 1954 1955 1956
## Darwin  1.4  0.9  0.6 -0.8  0.2 -1.6
## Tahiti  0.1 -1.1 -0.1  1.5  1.8  0.1

#df = data.frame(ptada)
dim(ptada)

## [1]  2 65

svdptd = svd(ptada) #SVD for the 2-by-65 matrix
U=round(svdptd$u, digits=2)
U

##       [,1] [,2]
## [1,] -0.75 0.66
## [2,]  0.66 0.75

D=round(diag(svdptd$d), digits=2)
D

##       [,1] [,2]
## [1,] 10.02 0.00
## [2,]  0.00 7.09

V =round(svdptd$v, digits=2)
t(V)[,1:5] #The first five year from 1951-1955

##       [,1]  [,2]  [,3] [,4] [,5]
## [1,] -0.10 -0.14 -0.05 0.16 0.10
## [2,]  0.14 -0.03  0.05 0.08 0.21

plot(1951:2015, -V[,1], type = 'o', 
     ylab = "PC", xlab = 'Year',
     main = "Tahiti-Darwin SLP Principal Components",
     col = 'black', ylim = c(-0.45, 0.45))
lines(1951:2015, -V[,2], type = 'l', 
      lty = 2, col = 'purple')
legend(1950, 0.5, lwd = 2, c("PC1", "PC2"), col = c('blue', 'Purple'),
       lty = c(1,2), bty="n")
#El Nino samples
points(c(1982, 1997), c(-V[32,1], -V[47,1]), 
       pch = 16, col = 'red')
text(1982, -V[32,1] + 0.07, "EN", col = 'red')
text(1997, -V[47,1] + 0.07, "EN", col = 'red')
#La Nina samples
points(c(1975, 2010), c(-V[25,1], -V[60,1]), 
       pch = 16, col = 'blue')
text(1975, -V[25,1] - 0.07, "LN", col = 'blue')
text(2010, -V[60,1] - 0.07, "EN", col = 'blue')

t = 1951:2015
y = -3*V[,1] + 0.3*cumsum(V[,1]) - 0.3*cumsum(V[,2])
plot(t, y, type = 'o', 
     ylab = "PC", xlab = 'Year',
     main = "Tahiti-Darwin SLP Principal Components",
     col = 'black', ylim = c(-1, 1))
lines(1951:2015, -V[,2], type = 'l', 
      lty = 2, col = 'purple')
legend(1950, 0.5, lwd = 2, c("PC1", "PC2"), col = c('blue', 'Purple'),
       lty = c(1,2), bty="n")
#El Nino samples
points(c(1982, 1997), c(-V[32,1], -V[47,1]), 
       pch = 16, col = 'red')
text(1982, -V[32,1] + 0.07, "EN", col = 'red')
text(1997, -V[47,1] + 0.07, "EN", col = 'red')
#La Nina samples
points(c(1975, 2010), c(-V[25,1], -V[60,1]), 
       pch = 16, col = 'blue')
text(1975, -V[25,1] - 0.07, "LN", col = 'blue')
text(2010, -V[60,1] - 0.07, "EN", col = 'blue')

t= 1951:2015
plot(t, cumsum(V[,1]), type = 'l', 
     ylim = c(-1,1), ylab = "Cumsum index")
lines(t, cumsum(V[,2]), 
      type = 'l', col = 'red') 
lines(t, -cumsum(V[,1]) + cumsum(V[,2]), 
      type = 'l', col = 'blue') 

#For a better computer display
library(plotly)

## Warning: package 'plotly' was built under R version 4.4.2

## Loading required package: ggplot2

## Warning: package 'ggplot2' was built under R version 4.4.2

## 
## Attaching package: 'ggplot2'

## The following objects are masked from 'package:psych':
## 
##     %+%, alpha

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

PC = data.frame(1951:2015, V)
colnames(PC) = c("Year", "PC1", "PC2")
plot_ly(data = PC, x = ~Year, y = ~(-PC1), 
        mode = 'lines+markers') %>%
  layout(yaxis = list(title = 'PC1'))

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter

Minors and Co-factor

minors <- function(b){
  n <- nrow(b)
  a <- matrix(NA, n, n)
  for(i in 1:n)
    for(j in 1:n)
      a[i, j] = det(b[-i, -j])
  a
}
b = matrix(c(2,3,5,6,7,1,9,4,5), 
           nrow = 3, ncol = 3)
minors(b) 

##      [,1] [,2] [,3]
## [1,]   31   -5  -32
## [2,]   21  -35  -28
## [3,]  -39  -19   -4

cofactors <- function(b) (-1)^(row(b)+col(b)) *minors(b)
cofactors(b)

##      [,1] [,2] [,3]
## [1,]   31    5  -32
## [2,]  -21  -35   28
## [3,]  -39   19   -4

#Adjoint matrix, aka, adjugate of a square matrix
A <- matrix(c(1,4,5,3,7,2,2,8,3),nrow=3,ncol=3)
A

##      [,1] [,2] [,3]
## [1,]    1    3    2
## [2,]    4    7    8
## [3,]    5    2    3

#install.packages('RConics')
library(RConics)

## 
## Attaching package: 'RConics'

## The following object is masked from 'package:psych':
## 
##     polar

## The following object is masked from 'package:pracma':
## 
##     polar

B <- adjoint(A)
B

##      [,1] [,2] [,3]
## [1,]    5   -5   10
## [2,]   28   -7    0
## [3,]  -27   13   -5

#B = det(A) * inverse of A
C = det(A)*solve(A)
C - B

##               [,1]         [,2]          [,3]
## [1,]  8.881784e-16 0.000000e+00  3.552714e-15
## [2,]  1.065814e-14 0.000000e+00 -7.195890e-16
## [3,] -7.105427e-15 1.776357e-15 -1.776357e-15

Characteristic polynomials, Cayley-Hamilton theorem, Jordon forms

A = matrix(c(1,1,0,2), nrow = 2)
A

##      [,1] [,2]
## [1,]    1    0
## [2,]    1    2

det(A)

## [1] 2

library(pracma)
cpA = charpoly(A, info = TRUE)

## Error term: 4

cpA

## $cp
## [1]  1 -3  2
## 
## $det
## [1] 2
## 
## $inv
##      [,1] [,2]
## [1,]  1.0  0.0
## [2,] -0.5  0.5

cpA$cp #1 x^2 - 3 x + 2

## [1]  1 -3  2

cpA$inv #= solve(A)

##      [,1] [,2]
## [1,]  1.0  0.0
## [2,] -0.5  0.5

cpA$det #= det(A)

## [1] 2

cpA$inv %*% A #zapsmall(cpA$inv %*% A), is equal to a 2-by-2 identity matrix

##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

#Cayley-Hamilton theorem
#cp(x) = a2 x^2 + a1 x + a0
#cp(A) = a2 A^2 + a1 A + a0 I = 0
library(expm)

## Warning: package 'expm' was built under R version 4.4.2

## 
## Attaching package: 'expm'

## The following object is masked from 'package:Matrix':
## 
##     expm

## The following objects are masked from 'package:pracma':
## 
##     expm, logm, sqrtm

A%^%2 - 3*A%^%1 + 2*diag(2) 

##      [,1] [,2]
## [1,]    0    0
## [2,]    0    0

#Jordan normal form
#install.packages("pracma")
library(pracma)
A = matrix(c(2, 1, 0, 2), nrow = 2)
A

##      [,1] [,2]
## [1,]    2    0
## [2,]    1    2

#jordan(A)

#install.packages("pracma")
library(pracma)
A <- matrix(c(2, 2, 0, 2), nrow = 2, byrow = TRUE)
#J <- jordan(A)
#J


#Jordan matrix
#install.packages('mcompanion')
library(mcompanion)

## Warning: package 'mcompanion' was built under R version 4.4.2

Jordan_matrix(4, 2)

##      [,1] [,2]
## [1,]    4    1
## [2,]    0    4

eigen(Jordan_matrix(4, 2))

## eigen() decomposition
## $values
## [1] 4 4
## 
## $vectors
##      [,1]          [,2]
## [1,]    1 -1.000000e+00
## [2,]    0  8.881784e-16

#$values: repeated eigenvalues 2 times

Jordan_matrix(5, 3)

##      [,1] [,2] [,3]
## [1,]    5    1    0
## [2,]    0    5    1
## [3,]    0    0    5

eigen(Jordan_matrix(5, 3))

## eigen() decomposition
## $values
## [1] 5 5 5
## 
## $vectors
##      [,1]          [,2]          [,3]
## [1,]    1 -1.000000e+00  1.000000e+00
## [2,]    0  1.110223e-15 -1.110223e-15
## [3,]    0  0.000000e+00  1.232595e-30

#$values: repeated eigenvalues 3 times

Jordan_matrix(6, 1)

##      [,1]
## [1,]    6

Jordan_matrix(6, 4)

##      [,1] [,2] [,3] [,4]
## [1,]    6    1    0    0
## [2,]    0    6    1    0
## [3,]    0    0    6    1
## [4,]    0    0    0    6

eigen(Jordan_matrix(6, 4))

## eigen() decomposition
## $values
## [1] 6 6 6 6
## 
## $vectors
##      [,1]          [,2]          [,3]          [,4]
## [1,]    1 -1.000000e+00  1.000000e+00 -1.000000e+00
## [2,]    0  1.332268e-15 -1.332268e-15  1.332268e-15
## [3,]    0  0.000000e+00  1.774937e-30 -1.774937e-30
## [4,]    0  0.000000e+00  0.000000e+00  2.364691e-45

#$values: repeated eigenvalues 4 times

A = Jordan_matrix(4, 3)
A

##      [,1] [,2] [,3]
## [1,]    4    1    0
## [2,]    0    4    1
## [3,]    0    0    4

cpA = charpoly(A, info = TRUE)

## Error term: 0

cpA #$cp

## $cp
## [1]   1 -12  48 -64
## 
## $det
## [1] 64
## 
## $inv
##      [,1]    [,2]      [,3]
## [1,] 0.25 -0.0625  0.015625
## [2,] 0.00  0.2500 -0.062500
## [3,] 0.00  0.0000  0.250000

library(expm)
A%^%3 - 12*A%^%2 + 48*A%^%1- 64*diag(3) 

##      [,1] [,2] [,3]
## [1,]    0    0    0
## [2,]    0    0    0
## [3,]    0    0    0

#Cayley-Hamilton theorem

## a matrix with the above 3 blocks
Jordan_matrix(c(4, 5, 6), c(2, 3, 1))

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    4    1    0    0    0    0
## [2,]    0    4    0    0    0    0
## [3,]    0    0    5    1    0    0
## [4,]    0    0    0    5    1    0
## [5,]    0    0    0    0    5    0
## [6,]    0    0    0    0    0    6

## a matrix with a 2x2 Jordan block for eval 1 and two simple 0 eval's
m <- make_mcmatrix(eigval = c(1), co = cbind(c(1,1,1,1), c(0,1,0,0)),
                   dim = 4, len.block = c(2))
m

##      [,1] [,2] [,3] [,4]
## [1,]    0    1    0    0
## [2,]   -1    2    0    0
## [3,]    0    1    0    0
## [4,]    0    1    0    0

m.X <- cbind(c(1,1,1,1), c(0,1,0,0), c(0,0,1,0), c(0,0,0,1))
m.X

##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    1    1    0    0
## [3,]    1    0    1    0
## [4,]    1    0    0    1

m.J <- cbind(c(1,0,0,0), c(1,1,0,0), rep(0,4), rep(0,4))
m.J

##      [,1] [,2] [,3] [,4]
## [1,]    1    1    0    0
## [2,]    0    1    0    0
## [3,]    0    0    0    0
## [4,]    0    0    0    0

from_Jordan(m.X, m.J) # == m

##      [,1] [,2] [,3] [,4]
## [1,]    0    1    0    0
## [2,]   -1    2    0    0
## [3,]    0    1    0    0
## [4,]    0    1    0    0

m.X %*% m.J %*% solve(m.X) # == m

##      [,1] [,2] [,3] [,4]
## [1,]    0    1    0    0
## [2,]   -1    2    0    0
## [3,]    0    1    0    0
## [4,]    0    1    0    0

all(m == from_Jordan(m.X, m.J)) && all(m == m.X %*% m.J %*% solve(m.X))

## [1] TRUE

## TRUE
## which column(s) in m.X correspond to 1st Jordan block?
chain_ind(1, c(2,1,1)) # c(1, 2) since 2x2 Jordan block

## [1] 1 2

## which column(s) in m.X correspond to 2nd Jordan block?
chain_ind(2, c(2,1,1)) # 3, simple eval

## [1] 3

## which column(s) in m.X correspond to 1st and 2nd Jordan blocks?
chain_ind(c(1, 2), c(2,1,1)) # c(1,2,3)

## [1] 1 2 3

## non-contiguous subset are ok:
chain_ind(c(1, 3), c(2,1,1)) # c(1,2,4)

## [1] 1 2 4

## split the chains into a list of matrices
chains_to_list(m.X, c(2,1,1))

## [[1]]
##      [,1] [,2]
## [1,]    1    0
## [2,]    1    1
## [3,]    1    0
## [4,]    1    0
## 
## [[2]]
##      [,1]
## [1,]    0
## [2,]    0
## [3,]    1
## [4,]    0
## 
## [[3]]
##      [,1]
## [1,]    0
## [2,]    0
## [3,]    0
## [4,]    1

m.X %*% m.J

##      [,1] [,2] [,3] [,4]
## [1,]    1    1    0    0
## [2,]    1    2    0    0
## [3,]    1    1    0    0
## [4,]    1    1    0    0

m %*% m.X # same

##      [,1] [,2] [,3] [,4]
## [1,]    1    1    0    0
## [2,]    1    2    0    0
## [3,]    1    1    0    0
## [4,]    1    1    0    0

all(m.X %*% m.J == m %*% m.X) # TRUE

## [1] TRUE

m %*% c(1,1,1,1) # = c(1,1,1,1), evec for eigenvalue 1

##      [,1]
## [1,]    1
## [2,]    1
## [3,]    1
## [4,]    1

m %*% c(0,1,0,0) # gen.e.v. for eigenvalue 1

##      [,1]
## [1,]    1
## [2,]    2
## [3,]    1
## [4,]    1

## indeed:
all( m %*% c(0,1,0,0) == c(0,1,0,0) + c(1,1,1,1) ) # TRUE

## [1] TRUE

## m X = X jordan.block
cbind(c(1,1,1,1), c(0,1,0,0)) %*% cbind(c(1,0), c(1,1))

##      [,1] [,2]
## [1,]    1    1
## [2,]    1    2
## [3,]    1    1
## [4,]    1    1

m %*% cbind(c(1,1,1,1), c(0,1,0,0))

##      [,1] [,2]
## [1,]    1    1
## [2,]    1    2
## [3,]    1    1
## [4,]    1    1