行列分解（関係データ学習　講談社、画像処理の統計モデリング　共立出版）

```{r}

#行列分解

#1次勾配降下法　p.101

library(dplyr)

n=100;m=10

X=matrix(sample(1:10,n*m,replace = T),nrow=n,ncol=m)

r=8

U=matrix(sample(1,n*r,replace = T),nrow=n,ncol=r)

V=matrix(sample(1,m*r,replace=T),nrow=m,ncol=r)

ite=500000

eta=0.0001

lambda=1

f_l2=3000000

option=0

cost=c()

for(j in 1:ite){
  
f_l2_pre=f_l2  
  
f_l2=sum((X-U%*%t(V))^2)/2+(lambda/2)*(sum(U^2)+sum(V^2)) 

if(option==1){
  
du=-(X-U%*%t(V))%*%V+lambda*U 
  
du=-eta*du

dv=-t(X)%*%U+V%*%(t(U)%*%U+lambda*diag(1,r))

dv=-eta*dv

U=U+du

V=V+dv

}else{
  
U=(1-eta)*U+eta*X%*%V%*%solve(t(V)%*%V+lambda*diag(1,r))  
  
V=(1-eta)*V+eta*t(X)%*%U%*%solve(t(U)%*%U+lambda*diag(1,r))  
  
}

#print(sum(abs(X-U%*%t(V))))

cost=c(cost,sum(abs(X-U%*%t(V))))

}

print(sum(abs(X-U%*%t(V))))


```
```{r}

#1次勾配降下法　p.101

#欠損値穴埋め

library(dplyr)
n=100;m=10

nums=17

n1=sample(1:n,nums,replace = T)
n2=sample(1:m,nums,replace = T)

M=array(1,dim=c(n,m))

for(i in 1:length(n1)){

M[n1[i],n2[i]]=0

}

X=matrix(sample(1:10,n*m,replace = T),nrow=n,ncol=m)
r=8

U=matrix(sample(1,n*r,replace = T),nrow=n,ncol=r)
V=matrix(sample(1,m*r,replace=T),nrow=m,ncol=r)

ite=500000

eta=0.001

lambda=1

f_l2=3000000

option=0

cost=c()

for(j in 1:ite){

f_l2_pre=f_l2  

f_l2=sum((X*M-U%*%t(V)*M)^2)/2+(lambda/2)*(sum(U^2)+sum(V^2)) 

if(option==1){

du=-(X*M-U%*%t(V)*M)%*%V+lambda*U 

du=-eta*du

dv=-t(X*M)%*%U+V%*%(t(U)%*%U+lambda*diag(1,r))

dv=-eta*dv

U=U+du

V=V+dv

}else{

U=(1-eta)*U+(eta*M*X)%*%V%*%solve(t(V)%*%V+lambda*diag(1,r))  

V=(1-eta)*V+eta*t(X*M)%*%U%*%solve(t(U)%*%U+lambda*diag(1,r))  

}

#print(sum(abs(X-U%*%t(V))))

cost=c(cost,sum(abs(X*M-U%*%t(V)*M)))

}

print(sum(abs(X*M-U%*%t(V)*M)))

repair=0

if(repair==1){

for(i in 1:length(n1)){
  
X[n1[i],n2[i]]=sum(U[n1[i],]*V[n2[i],])  
  
}

}  
  
print(sum(abs(X-U%*%t(V))))

```

```{r}

#CP分解

data=data.frame(i=0,j=0,k=0,value=0)

I=10;J=12;K=20;r=9

eta=0.001

for(i in 1:I){
for(j in 1:J){  
for(k in 1:K){ 
 
data_sub=data.frame(i=i,j=j,k=k,value=rpois(1,sample(10:100,1)))    
  
data=rbind(data,data_sub)

}
}
}

data=tail(data,nrow(data)-1)

U=array(1,dim=c(I,r))

V=array(2,dim=c(J,r))

W=array(3,dim=c(K,r))


#mat1=V[,r]%*%t(W[,r])

#mat2=U[,r]%*%t(W[,r])

#mat3=U[,r]%*%t(V[,r])

sita1=array(1,dim=c(I,r))

sita2=array(1,dim=c(J,r))

sita3=array(1,dim=c(K,r))

for(i in 1:I){
  
X1=array(0,dim=dim(V[,r]%*%t(W[,r])))  
  
for(s in 1:r){
  
mat1=V[,s]%*%t(W[,s])

X1=X1+U[i,s]*mat1  
  
}

for(s in 1:r){

sita1[i,s]=t(V[,s])%*%X1%*%W[,s]

}


}


for(j in 1:J){
  
X2=array(0,dim=dim(U[,r]%*%t(W[,r])))  
  
for(s in 1:r){
  
mat2=U[,s]%*%t(W[,s])

X2=X2+V[j,s]*mat2  
  
}

for(s in 1:r){

sita2[j,s]=t(U[,s])%*%X2%*%W[,s]

}


}


for(k in 1:K){
  
X3=array(0,dim=dim(U[,r]%*%t(V[,r])))  
  
for(s in 1:r){
  
mat3=U[,s]%*%t(V[,s])

X3=X3+W[k,s]*mat3 
  
}

for(s in 1:r){

sita3[k,s]=t(U[,s])%*%X3%*%V[,s]

}


}

for(i in 1:I){
  
U[i,]=U[i,]-  
  
}

```

```{r}

#Healy-Westmacott Procedure as an EM algorithm p.55

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

X=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

Y=data$y

n=nrow(X);m=18

Y=Y[1:m]

beta=rep(1,ncol(X))

sigma=1

ite=100

for(l in 1:ite){

y_sub=X[(m+1):n,]%*%beta

beta=solve(t(X)%*%X)%*%t(X)%*%c(Y,y_sub)

sigma=(sum((Y-X[1:m,]%*%beta)^2)+(n-m)*sigma)/n

print(sigma)

}

```

```{r}

#Buck`s method p.51

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

X=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

Y=data$y

n=nrow(X)

nums=20;W=X;

n1=sample(1:dim(X)[1],nums,replace=T)
n2=sample(1:dim(X)[2],nums,replace=T)

W[n1,n2]=NA;W_sub=W

beta=c(1,1,1)

W_ave=c()

for(j in 1:ncol(W)){
  
vec=W[,j];  
  
W_ave=c(W_ave,mean(vec[is.na(vec)!=T]))  
  
}

ite=100


for(s in 1:ite){

W_pre=W  
    
W=W_sub  

for(j in 1:nrow(W)){
for(i in 1:ncol(W)){
 
  
if(is.na(W[j,i])>0){  
  
W_vec=(t(W)-c(W_ave));W_vec[is.na(W_vec)==T]=0

W_vec=apply(W_vec,1,sum)  
  
W[j,i]=W_ave[i]+sum(beta*W_vec)  
  
}  
  
}  
}
  
beta=solve(t(W)%*%W)%*%t(W)%*%Y  

print(sum(abs(W_pre-W)))

W_ave=apply(W,2,mean)
  
}


```

```{r}

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

X=as.matrix(data[,colnames(data) %in% c("x1")])

Y=as.numeric(data$y)

f1=function(x){
  
return((x^2))  
  
}

f2=function(x){
  
return(x)  
  
}

f_mat=cbind(f1(X),f2(X))



eta=0.01

a=array(1,dim=c(2,length(Y)))

ite=10000

for(i in 1:ite){
  
a_pre=a  
  
H=a/(t(f_mat)%*%f_mat%*%a)
  
a=a+H*(t(f_mat)%*%(Y-f_mat%*%a))

print(sum(abs(Y-f_mat%*%a)))
  
}




```


```{r}

#f fixed

library(dplyr)

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

X=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

Y=matrix(sample(data$y,nrow(data)^2,replace=T),ncol=nrow(data))



f_mat=X

eta=0.01

a=array(1,dim=c(ncol(X),nrow(data)))

ite=10000

for(i in 1:ite){
  
a_pre=a  
  
H=a/(t(f_mat)%*%f_mat%*%a)
  
a=a+H*(t(f_mat)%*%(Y-(f_mat%*%a)))

print(sum(abs(Y-f_mat%*%a)))
  
}




```

```{r}

#A fixed

library(dplyr)


data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22))

X=as.matrix(data[,colnames(data) %in% c("x1","x2","x3")])

Y=matrix(sample(data$y,nrow(data)^2,replace = T),ncol=nrow(data))



f_mat=X

eta=0.01

a=array(sample(1:10,ncol(X)*nrow(X),replace=T),dim=c(ncol(X),nrow(X)))

ite=1000

for(i in 1:ite){
  
a_pre=a  
  
H=f_mat/(t(f_mat%*%a)%*%t(a))

#H=f_mat/(f_mat%*%a%*%t(a))
  
f_mat=f_mat+H*(t(Y-f_mat%*%a)%*%t(a))

print(sum(abs(Y-f_mat%*%a)))
  
}




```

```{r}

#因子分析　-その理論と方法-

#因子分析における推定(回帰と行列分解)　p.77

M=20

n=100;m=20

X=array(sample(1:20,n*m,replace = T),dim=c(n,m))

Y=sample(20:40,n,replace=T)

U=diag(1,m)

S=cov(X)

eigen_values=eigen(solve(U)%*%S%*%solve(U))$values[1:M]

eigen_vectors=eigen(solve(U)%*%S%*%solve(U))$vectors[,1:M]

for(j in 1:ncol(eigen_vectors)){

eigen_vectors[,j]=eigen_vectors[,j]/sum(eigen_vectors[,j]^2)

}

omega=eigen_vectors

lambda=U%*%omega%*%diag(sqrt(eigen_values))

f=X%*%solve(U)%*%omega%*%diag(1/sqrt(eigen_values))

t=f%*%t(lambda)

tau2=sum(diag((X-t)%*%solve(U^2)%*%t(X-t)))/(n*m)

sum(abs(X-t))

```