(工業基礎)最適化とその応用　矢部博先生　数理工学社

```r

#p.108
#主双対内点法

A=t(matrix(c(2,1,-3,5,3,-7,4,1),nrow=4))
C=c(3,2,5,1)
b=c(10,15)

x=rep(1,4)
z=rep(1,4)
y=rep(1,2)
gamma=1 #近傍系のガンマ。
mesh=0.001 #学習率の刻み
mu=sum(x*z)/length(z)
tau=mu*gamma;tau=tau-tau/10
sig=tau/mu
learning_max=-sqrt(2)*(sig-gamma)/((1+1/gamma)*length(x))
learning_rate_vec=seq(mesh,learning_max,mesh) #学習率に幅を持たせる
ite=1000 #反復計算回数
Neighbourhood=c()
Neighbourhood_No=c()

for(l in 1:ite){

x_pre=x
y_pre=y
z_pre=z
  
mu=sum(x*z)/length(z)
sig=tau/mu

mat1=cbind(diag(0,length(x)),t(A),diag(1,length(z)))
mat2=cbind(A,diag(0,length(y)),array(0,dim=c(length(y),length(z))))
mat3=cbind(diag(z),array(0,dim=dim(t(A))),diag(x))
mat=rbind(mat1,mat2,mat3)

r=-c(t(A)%*%y+z-C,A%*%x-b,diag(x)%*%diag(z)%*%rep(1,length(x))-sig*mu*rep(1,length(x)))

dxyz=solve(mat)%*%r
xyz_mat=dxyz%*%learning_rate_vec+c(x,y,z)
xyzsum=apply(ifelse(xyz_mat>0,1,0)[-c((length(x)+1):(length(x)+2)),],2,sum)
norm_dxdz=sqrt(sum((dxyz[c(1:length(x))]*dxyz[((length(c(x,y)))+1):length(c(x,y,z))])^2))


Neighbourhood_bounded=length(x)*mu*(1+1/gamma)/(2*sqrt(2))
select_xyz_column=c(1:length(learning_rate_vec))[xyzsum==length(c(x,z))]
xyz=xyz_mat[,max(select_xyz_column)]

x=xyz[1:length(x)];xyz=xyz[-c(1:length(x))]
y=xyz[1:length(y)];xyz=xyz[-c(1:length(y))]
z=xyz

if(norm_dxdz<=Neighbourhood_bounded){
Neighbourhood=c(Neighbourhood,xyz_mat[,max(select_xyz_column)]) 
Neighbourhood_No=c(Neighbourhood_No,l)
}

if(tau>10^(-10)){
tau=tau/log(l+1) #反復計算により、小さくしていく
}

print(mu)

}

#近傍系(収束列が取れる)
Neighbourhood=data.frame(matrix(Neighbourhood,nrow=length(c(x,y,z))))
colnames(Neighbourhood)=Neighbourhood_No

#主問題の実行可能性
A%*%x-b

#双対問題の実行可能性
t(A)%*%y+z-C

#相補性条件
sum(z*x)

```

```r

#共役勾配法 p.146

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

f=function(x){
x=c(x)
return(as.numeric(x%*%A%*%x)/2-sum(c(6,7)*x))
}

df=function(x){
dfval=c()
for(j in 1:length(x)){
xsub=x;xsub[j]=xsub[j]+h  
dfval=c(dfval,(f(xsub)-f(x))/h)  
}
return(dfval)
}

Armijo_condition=1 #アルミホルールを使用する

#betaの指定方法
fletcher_reeves=0
polak_ribiere=0
hestenes_stiefel=1
dai_yuan=0


guzai=10^(-4)
x=c(2,1)
h=10^(-3)
l=1
x_pre=x+1
beta=1
beta_vec=c(1,exp(cumsum(rep(log(tau),1000))))
alpha=0.1

while(sum(abs(f(x_pre)-f(x)))>10^(-9) & alpha>0)({

f_pre=f(x)
x_pre=x
beta_pre=beta
if(l==1){  
d=-df(x)  
}else{
d=-df(x)+beta*d 
}
 

if(Armijo_condition==1){
assign_no=c()
for(j in 1:length(beta_vec)){
if(f(x+beta_vec[j]*d)<=f(x)+guzai*beta_vec[j]*sum(df(x)*d)){
assign_no=c(assign_no,j)  
}}
if(length(assign_no)==0){
alpha=0  
}else{
alpha=max(beta_vec[assign_no])
}
}else{
alpha=sum(d*(-df(x)))/as.numeric(d%*%A%*%d) 
}

df_pre=df(x) 
x=x+alpha*d
y=df(x)-df_pre

if(fletcher_reeves==1){
beta=sum(df(x)^2)/sum(df_pre^2)
}

if(polak_ribiere==1){
beta=sum(df(x)*y)/sum(df_pre^2)  
}

if(hestenes_stiefel==1){
beta=sum(df(x)*y)/sum(d*y)  
}

if(dai_yuan==1){
beta=sum(df(x)^2)/sum(d*y)  
}

print(f(x))
l=l+1

})


```

```r

#準ニュートン法 p.156

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

f=function(x){
x=c(x)
return(as.numeric(x%*%A%*%x)/2-sum(c(6,7)*x))
}

df=function(x){
dfval=c()
for(j in 1:length(x)){
xsub=x;xsub[j]=xsub[j]+h  
dfval=c(dfval,(f(xsub)-f(x))/h)  
}
return(dfval)
}

#Bの更新
DFP=0　#Davidon-Fletcher-Powellの公式
BFGS=1　#Broyden-Fletcher-Goldfarb-Shannoの公式
SROU=0 #symmetric rank one update(対象ランクワン公式)
SM=0 #Sherman-Morrisonの公式(BFGSの逆行列)

tau=0.94
guzai=10^(-4)
x=c(2,1)
h=10^(-3)
l=1
beta=1
beta_vec=c(1,exp(cumsum(rep(log(tau),1000))))
alpha=0.1
B=diag(1,length(x))
x_pre=x+1
d=-df(x)

while(sum(abs(f(x_pre)-f(x)))>10^(-9) & alpha>0)({

f_pre=f(x)
x_pre=x
beta_pre=beta

#Armijoのルール
assign_no=c()
for(j in 1:length(beta_vec)){
if(f(x+beta_vec[j]*d)<=f(x)+guzai*beta_vec[j]*sum(df(x)*d)){
assign_no=c(assign_no,j)  
}}
if(length(assign_no)==0){
alpha=0  
}else{
alpha=max(beta_vec[assign_no])
}

d=-solve(B)%*%df(x)
df_pre=df(x) 
x=x+alpha*d
y=df(x)-df_pre
s=x-x_pre

if(BFGS==1){
B=B-as.numeric(B%*%s%*%t(B%*%s))/as.numeric(c(s)%*%B%*%s)+y%*%t(y)/sum(s*y)
}

if(DFP==1){
B=B-((B%*%s)%*%t(y)+y%*%t(B%*%s))/sum(s*y)+(1+as.numeric(c(s)%*%B%*%c(s))/sum(s*y))*y%*%t(y)/sum(s*y)
}

if(SROU==1){
B=B+(y-B%*%s)%*%t(y-B%*%s)/sum(s*c(y-B%*%s))  
}

if(SM==1){
H=solve(B)
H=H-(H%*%y%*%t(s)+s%*%t(H%*%y))/sum(s*y)+(1+as.numeric(c(y)%*%H%*%c(y))/sum(s*y))*s%*%t(s)/sum(s*y)
B=solve(H)
}

print(f(x))
l=l+1

})




```

```r

#ガンマ分布の高速最適化アルゴリズム
#局所近似最大化法

#サンプルデータ
N = 500
values = rgamma(N, 3, scale = 2)

#初期パラメータと反復計算回数の指定
a=5 #初期値
b=5 #初期値
ite=10　#反復計算回数

for(l in 1:ite){
mv = mean(values)
a = 1/ (1/a + (mean(log(values)) - digamma(a) - log(mv) + log(a)) / (a^2 * (1/a - trigamma(a))))
b = mv / a
log_lik = sum(log(dgamma(values, a, b)))
print(paste0("a:",round(a,3),",b:",round(b,3)))
}  

```

```r


#準ニュートン法での最小化 
#ガンマ分布の最適化アルゴリズム

#ガンマ分布の場合
サンプルデータ
N = 500
values = rgamma(N, 3, scale = 2)

f=function(x){　#最小化したい関数
a=x[1];b=x[2]
logliklihood=(a-1)*sum(log(values))-N*log(gamma(a))-N*a*log(b)-sum(values)/b
return(-logliklihood)
}

df=function(x){　#関数の勾配ベクトル
dfval=c()
for(j in 1:length(x)){
xsub=x;xsub[j]=xsub[j]+h  
dfval=c(dfval,(f(xsub)-f(x))/h)  
}
return(dfval)
}

#Bの更新
DFP=0　#Davidon-Fletcher-Powellの公式
BFGS=0　#Broyden-Fletcher-Goldfarb-Shannoの公式
SROU=1 #symmetric rank one update(対象ランクワン公式)
SM=0 #Sherman-Morrisonの公式(BFGSの逆行列)

x=c(2,1) #初期値
h=10^(-3)　#差分幅
B=diag(1,length(x)) #近似行列(正定値対称)
x_pre=x+1　#初期値
d=-df(x) 
alpha=0.0008 #学習幅(ガンマ分布の場合)
ite=10000　#反復回数(ガンマ分布の場合)

for(s in 1:ite){

f_pre=f(x)
x_pre=x

d=-solve(B)%*%df(x)
df_pre=df(x) 
x=x+alpha*d
y=df(x)-df_pre
s=x-x_pre

if(BFGS==1){
B=B-as.numeric(B%*%s%*%t(B%*%s))/as.numeric(c(s)%*%B%*%s)+y%*%t(y)/sum(s*y)
}

if(DFP==1){
B=B-((B%*%s)%*%t(y)+y%*%t(B%*%s))/sum(s*y)+(1+as.numeric(c(s)%*%B%*%c(s))/sum(s*y))*y%*%t(y)/sum(s*y)
}

if(SROU==1){
B=B+(y-B%*%s)%*%t(y-B%*%s)/sum(s*c(y-B%*%s))  
}

if(SM==1){
H=solve(B)
H=H-(H%*%y%*%t(s)+s%*%t(H%*%y))/sum(s*y)+(1+as.numeric(c(y)%*%H%*%c(y))/sum(s*y))*s%*%t(s)/sum(s*y)
B=solve(H)
}

print(f(x))　#最小化したい関数の値

}

#反復式での解
x

```


```r

#p.224 主双対内点法(確認中)

data("EuStockMarkets")
mat=as.matrix(EuStockMarkets)[-1,]/1000


initial_rate=EuStockMarkets[1,]
C=10000

g=function(x){
return((sum(x*initial_rate)-C)^2)  
}

f=function(x){
return(-sum(mat%*%x))  
}

L=function(x,mu,lam){
return(f(x)-mu*sum(log(x))+lam*g(x))  
}

h=0.001
Mc=10^5
tau=1-10^(-5)
alpha=0.01

param=rep(1,4)
z=rep(1,4)
MU=10
LAM=5
ite=4*10^4

B=diag(1,length(param))

for(l in 1:ite){
  
if(l>1){
dL_pre=dL  
}

dL=c()
dg=c()
for(j in 1:length(param)){
param_sub=param
param_sub[j]=param_sub[j]+h
dL=c(dL,(L(param_sub,MU,LAM)-L(param,MU,LAM))/h)
dg=c(dg,(g(param_sub)-g(param))/h)
}

if(l==1){
q=dL  
s=param
}else{
q=dL-dL_pre
s=param-param_pre
}
q=t(t(c(q)))
s=c(s)

pthi=ifelse(as.numeric(s%*%q)>=0.2*as.numeric(s%*%B%*%s),1,0.8*as.numeric(s%*%B%*%s)/as.numeric(s%*%(B%*%s-q)))

B=B-(B%*%s)%*%t(B%*%s)/as.numeric(s%*%B%*%s)+q%*%t(q)/as.numeric(s%*%q)

J1=cbind(B,dg,-diag(1,length(param)))
J2=c(t(c(dg)),0,rep(0,length(param)))
J3=cbind(diag(c(z)),rep(0,length(param)),diag(c(param)))
J=rbind(J1,J2,J3)
r=c(dL,g(param),(param*z)-MU)

w=solve(J)%*%(-r)
param_pre=param
if(sum(c(c(param,LAM,z)+alpha*w)<0)>0){
alpha=-c(param,LAM,z)/w
alpha=min(alpha[w<0])
}
parameters=c(param,LAM,z)+alpha*w
parameters[parameters==0]=0.0001
param=parameters[1:4];parameters=parameters[-c(1:4)]
LAM=parameters[1]
z=parameters[-1]
MU=MU*tau
print(sqrt(sum(r^2)))

}

```