スパースモデリング　基礎から動的システムへの応用　永原正章先生　コロナ社

```{r}

#一般化lassoを解くためのADMMアルゴリズム p.65


data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),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))

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

X=cbind(rep(1,nrow(X)),X)

Y=data$y

m=5

pthi=matrix(sample(1:9,m*ncol(X),replace=T),ncol=ncol(X))/10

mu=1

lam=0.01

w=rep(1,ncol(X))

z=rep(1,m)

v=rep(1,m)

ite=10^4

cost_vec=c()

for(l in 1:ite){

w=solve(t(X)%*%X+(1/mu)*t(pthi)%*%pthi)%*%(t(X)%*%Y+(1/mu)*t(pthi)%*%(z-v))  
  
dw=pthi%*%w+v

z=ifelse(abs(dw)<=lam*mu,0,ifelse(dw>lam*mu,dw-lam*mu,dw+lam*mu))

v=v+pthi%*%w-z

cost=sum((Y-X%*%w)^2)+lam*sum(abs(pthi%*%w))

print(cost)

cost_vec=c(cost_vec,cost)

}

plot(c(1:length(cost_vec)))


```

```{r}


#p.79 l0最適化問題(4.1)を解くためのマッチング追跡アルゴリズム

n=5;m=5

X=matrix(sample(c(0:1),n*m,replace=T),ncol=n)

pthi=X

Y=sample(1:9,m,replace=T)

x=rep(0,ncol(pthi))

r=Y

ite=1000

for(l in 1:ite){
  
x_pre=x
  
value=(t(pthi)%*%r)^2/diag(t(pthi)%*%pthi)

ik=c(1:length(value))[value==max(value)]

ik=ifelse(length(ik)>1,ik[1],ik)
  
xk=sum(pthi[,ik]*r)/sum(pthi[,ik]^2)

x[ik]=x[ik]+xk  

r=r-xk*pthi[,ik]

print(sum(abs(Y-pthi%*%x)))
  
}



```

```{r}

#反復ハードしきい値アルゴリズム(IHT) p.84 と反復sスパースアルゴリズム p.87


data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),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))

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

X=cbind(rep(1,nrow(X)),X)

Y=data$y

#定理4.4にのっとって設定

r=0.9/eigen(t(X)%*%X)$values[1]

#閾値を設定する場合

s=1

lam=1

x=rep(1,ncol(X))

ite=10^5

for(l in 1:ite){
  
dx=x-r*t(X)%*%(X%*%x-Y)

if(s>0){
  
dx=ifelse(abs(dx)<s,0,dx)
  
}else{  

dx=ifelse(abs(dx)<sqrt(2*r*lam),0,dx)

}

x=dx

print(sum(Y-X%*%x))
  
}


```

```{r}

#L1最適化問題ADMMアルゴリズム　p.167

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

exps=function(s){
Ps=eigen(A*s)$vectors  
return((t(Ps)%*%diag(exp(eigen(A*s)$values))%*%Ps)%*%b)
}

t=1
n=5
h=t/n

Ah=A*h  
P=eigen(Ah)$vectors  
Ad=t(P)%*%diag(exp(eigen(Ah)$values))%*%P

#積分区間
aa=0;bb=h
#n;差分メッシュ
#積分の近似的計算：微分積分要論　学術図書出版　戸田暢茂先生　p.95の例4参照
N=10^1
mesh=aa+(bb-aa)*seq(1,N,1)/N
value=rep(0,length(b))
for(j in 1:N){
value=cbind(value,exps(mesh[j]))
}
#積分の結果 (p.82 (3.29))
bvec=apply(value[,-1],1,sum)*(bb-aa)/N

pthi=bvec
tueta=-Ad%*%x0
for(i in 1:(n-1)){
bvec=Ad%*%bvec  
pthi=cbind(bvec,pthi)
tueta=Ad%*%tueta
}

inv=solve(2*diag(1,n)+t(pthi)%*%(pthi))
M=inv%*%t(rbind(diag(1,n),diag(1,n),pthi))

#初期値
z=rep(1,2*n+length(b))
v=rep(1,2*n+length(b))/100

ite=100000
r=0.01
eta=0.01
for(l in 1:ite){
vpre=v;zpre=z
if(l>1){
u=u*(1-eta)+eta*M%*%(z-v)  
}else{
u=M%*%(z-v) 
}
v0=v[1:n]
v1=v[(n+1):(2*n)]
v2=v[-c(1:(2*n))]
z0=u+v0
z0[(u+v0)>=r]=(u+v0-r)[(u+v0)>=r]
z0[abs(u+v0)<r]=0
z0[(u+v0)<=(-r)]=(u+v0+r)[(u+v0)<=(-r)]
z1=sign(u+v1)*min(abs(u+v1),1)
z=z*(1-eta)+eta*c(z0,z1,tueta)
v=v*(1-eta)+eta*(rbind(diag(1,n),diag(1,n),pthi)%*%u-z)

print(sum(abs(z-zpre))+sum(abs(v-vpre)))

}

```