#library('matlib')
library('pracma')
#library('spatstat')

C29 pg. 354 Alt text

A = cbind(c(1,-3,1,1), c(9,-27,11,7), c(9,-29,13,7),  c(24,-68,26,18))
I = diag(nrow(A))
eig_vals = eig(A) 
A
##      [,1] [,2] [,3] [,4]
## [1,]    1    9    9   24
## [2,]   -3  -27  -29  -68
## [3,]    1   11   13   26
## [4,]    1    7    7   18
eig_vals
## [1] 2.000000e+00 2.000000e+00 1.000000e+00 3.064216e-14
# eigen values are 2, 2, 1, 0
# (A - lambdaI)x = 0, basis for eigenspace are eigenvectors, is null space
lambda = 2
rref_lambda_2 = rref(A - lambda*I)
rref_lambda_2
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0 -1.5
## [2,]    0    1    1  2.5
## [3,]    0    0    0  0.0
## [4,]    0    0    0  0.0
# eps_lambda_2 = N(A - 2*I)
# x3 and x4 are free variables
eig_v1 = c(0,-1, 1, 0) # when x3 is 1, other free variables are 0
eig_v2 = c(1.5, -2.5, 0, 1) # when x3 is 1, other free variables are 0
lambda = 1
rref_lambda_1 = rref(A - lambda*I)
rref_lambda_1
##      [,1] [,2] [,3]      [,4]
## [1,]    1    0    0 -1.666667
## [2,]    0    1    0  4.333333
## [3,]    0    0    1 -1.666667
## [4,]    0    0    0  0.000000
# eps_lambda_1 = N(A - 1*I)
# x4 is free variable
eig_v3 = c(5/3, -13/3, 5/3, 1) # set x4 to 1
lambda = 0
rref_lambda_0 = rref(A - lambda*I)
rref_lambda_0
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0   -3
## [2,]    0    1    0    5
## [3,]    0    0    1   -2
## [4,]    0    0    0    0
# eps_lambda_0 = N(A - 0*I)
# x4 is free variable
eig_v4 = c(3, -5, 2, 1)
# S = 4x4 matrix of eigvectors
S = cbind(eig_v1, eig_v2, eig_v3, eig_v4)
S
##      eig_v1 eig_v2    eig_v3 eig_v4
## [1,]      0    1.5  1.666667      3
## [2,]     -1   -2.5 -4.333333     -5
## [3,]      1    0.0  1.666667      2
## [4,]      0    1.0  1.000000      1
# A is diagonaizable iff geometric multiplicity == algebraic multiplicty for each eigenvalue
# characteristic polynomial:
# p_A(x) = (x-2)^2 * (x-1) * (x)
# when lambda = 2, geom multi == algeb multi == 2 (twice in polynomial; eig_v1 and eig_v2)
# when lambda = 1, geom multi == algeb multi == 1 (once in polynomial; eig_v3)
# when lambda = 0, gemo multi == algeb multi == 1 (once in polynomial; eig_v4)
# Thus A is diagonalizable 

# D = S-invAS
D = inv(S)%*%A%*%S
D
##               eig_v1        eig_v2        eig_v3        eig_v4
## eig_v1  2.000000e+00 -1.278977e-13 -1.278977e-13 -1.421085e-13
## eig_v2  0.000000e+00  2.000000e+00 -4.263256e-14 -2.842171e-14
## eig_v3 -2.842171e-14  4.973799e-14  1.000000e+00 -1.421085e-14
## eig_v4 -5.329071e-15 -4.329870e-15 -1.369275e-14 -1.576517e-14
# round to show nicer D
round(D)
##        eig_v1 eig_v2 eig_v3 eig_v4
## eig_v1      2      0      0      0
## eig_v2      0      2      0      0
## eig_v3      0      0      1      0
## eig_v4      0      0      0      0
# check SDS-inv = A
result = S%*%D%*%inv(S)
result
##      [,1] [,2] [,3] [,4]
## [1,]    1    9    9   24
## [2,]   -3  -27  -29  -68
## [3,]    1   11   13   26
## [4,]    1    7    7   18
A
##      [,1] [,2] [,3] [,4]
## [1,]    1    9    9   24
## [2,]   -3  -27  -29  -68
## [3,]    1   11   13   26
## [4,]    1    7    7   18