# http://javaspecialist.co.kr/board/795 게시글 참고

rm(list=ls())
# 다음 데이터 A를 특이값 분해 하세요.
# 설명참고 http://web.mit.edu/be.400/www/SVD/Singular_Value_Decomposition.htm
# A = U %*% S  %*% t(V)

(A <- matrix( c(2,1,0,0,4,3,0,0), ncol=2)) # 데이터

# 특이값과 U를 구하기 위해 A와 A의 전치행렬을 곱합니다.
sqrt(AAT <- A %*% t(A)) 
# 첫 번째 행이 -9.866069 14.000000    0    0
# 이것이 단위 벡터가 되게 하기 위해서는 

# R의 eigen() 함수는 고유값과 고유값 행렬을 구해줍니다.
# 고유값은 AAT - lambdaI = 0 을 이용해서도 구할 수 있습니다.
# det(AAT-LI) = 0 에서 L을 구하는데 그것이 람다값이다.
ev <- eigen(AAT)$values #29.8660687  0.1339313  0.0000000  0.0000000

W <- (A %*% t(A) - ev * diag(4))
W

eig_v1 <- c(1/sqrt(W[1,] %*% W[1, ])) * W[1,]
eig_v1 # 이것이 첫 번째 고유벡터

eig_v2 <- c(1/sqrt(W[1,] %*% W[1, ])) * W[2,]
eig_v2 # 이것이 두 번째 고유벡터

eig_v3 <- c(1/sqrt(W[1,] %*% W[1, ])) * W[3,]
eig_v3 # 이것이 세 번째 고유벡터

eig_v4 <- c(1/sqrt(W[1,] %*% W[1, ])) * W[4,]
eig_v4 # 이것이 네 번째 고유벡터

U <- matrix(c(eig_v2, eig_v1, eig_v3, eig_v4), ncol=4)
U

W <- (t(A) %*% A - ev * diag(2))
W
eig_v1 <- c(1/sqrt(W[1,] %*% W[1, ])) * W[1,]
eig_v1 # 이것이 첫 번째 고유벡터

eig_v2 <- c(1/sqrt(W[1,] %*% W[1, ])) * W[2,]
eig_v2 # 이것이 첫 번째 고유벡터

V <- matrix(c(eig_v2, eig_v1), ncol=2)
V


# AAT-LI에서 L이 0.134일 경우 x1, x2, ...를 구하면 
# 이것이 두 번째 고유벡터
# 특이값은 고유값에 루트를 씌워야 합니다.

(ev_sqrt <- sqrt(ev)) # 5.47, 0.37

# V를 구하기 위해 A의 전치행렬과 A행렬을 곱합니다.
ATA <- t(A) %*% A
eigen(ATA)$values #AAT를 이용해 구한 교유값과 동일한 값을 출력합니다.

# 특이값을 이용해 Sigma 행렬을 만들어줍니다.
(S <- diag(ev_sqrt))
(S <- S[,-c(3,4)])


# 이제 U, S, V를 얻었으니 증명을 해야 합니다.
# VT는 V의 전치행렬이며, AT는 A의 전치행렬, UT는 U의 전치행렬 입니다.
A
U %*% S %*% t(V)
#A=USVT and AT=VSUT
#ATA = VSUTUSVT
#ATA = VS2VT
#ATAV = VS2