非線形計画法 　オーム社  　矢部・八巻先生

```{r}
#unconstrained conlinear least squares problem

#Gauss-Newton method p.95

r=function(x){

x1=x[1];x2=x[2]

f1=x1^2-2*x2+1

f2=x1^2-x2^2

return(c(f1,f2))

}



X=sample(1:9,2)

h=0.001

ite=10^4

eta=0.001


for(l in 1:ite){

J=array(0,dim=c(length(X),length(X)))

for(j in 1:length(X)){

vec=X;vec[j]=vec[j]+h

J[,j]=(r(vec)-r(X))/h

}

J=t(J)

d=-solve(t(J)%*%J)%*%t(J)%*%r(X)

cost=sum((J%*%d+r(X))^2)

print(cost)

X=X+eta*d

}



```

```{r}

#4.2.3 構造化準ニュートン法を用いた解法

#p.97

r=function(x){

x1=x[1];x2=x[2]

f1=x1^2-2*x2+1

f2=x1^2-x2^2

return(c(f1,f2))

}

Jab=function(y){

J=array(0,dim=c(length(y),length(y)))

for(j in 1:length(y)){

vec=y;vec[j]=vec[j]+h

J[,j]=(r(vec)-r(y))/h

}

return(J)

}

X=sample(1:9,2)

h=0.001

ite=10^5

eta=10^(-3)

X_pre=rep(-1,2)

for(l in 1:ite){

J=t(Jab(X))

s=X-X_pre

q=t(Jab(X)-Jab(X_pre))%*%r(X)

A=diag(c(q/s))

d=-solve(t(J)%*%J+A)%*%t(J)%*%r(X)

X_pre=X

X=X+eta*d

print(r(X))

}



```

```{r}

#p.98 Dennis。サイジングの方法。(4.2.18)の式で計算。

f=function(x){

x1=x[1];x2=x[2]

f1=x1^2-x2

f2=x1-x2^2

return(c(f1,f2))

}

Jacob=function(y){

vec=y

J=array(0,dim=c(length(y),length(y)))

for(j in 1:length(y)){

vec_sub=vec;vec_sub[j]=vec_sub[j]+h

J[,j]=(f(vec_sub)-f(vec))/h

}

return(J)

}


x=rep(5,2)

x_pre=rep(0,2)

A=diag(rep(1,2))

ite=10^3

h=0.01

eta=10^(-2)

for(l in 1:ite){


s=x-x_pre

q=t(Jacob(x)-Jacob(x_pre))%*%f(x)

beta=min(c(abs(sum(s*q)/as.numeric(t(c(s))%*%A%*%t(t(c(s))))),1))

z=t(Jacob(x))%*%Jacob(x)%*%s+q

A=A*beta+(t(t(c(q-beta*A%*%s)))%*%t(c(z))+t(t(c(z)))%*%t(c(q-beta*A%*%s)))/sum(s*z)-sum(s*(q-beta*A%*%s))*t(t(c(z)))%*%t(c(z))/(sum(s*z)^2)

d=-solve(t(Jacob(x))%*%Jacob(x)+A)%*%t(Jacob(x))%*%f(x)

x_pre=x

x=x+eta*d

print(sum(abs(f(x))))

}

```

```{r}

#p.100 Yabe and Yamaki

f=function(s){

x1=s[1];x2=s[2]

f1=x1^2-x2

f2=x1-x2^2

return(c(f1,f2))

}

Jacob=function(s){

vec=s

J=array(0,dim=c(length(s),length(s)))

for(j in 1:length(s)){

vec_sub=vec;vec_sub[j]=vec_sub[j]+h

J[,j]=(f(vec_sub)-f(vec))/h

}

return(J)

}


ite=10^3

x=rep(5,2)

x_pre=rep(0,2)

eta=0.01

W=diag(rep(1,length(x)))

h=0.01

for(l in 1:ite){

s=x-x_pre

z=t(Jacob(x))%*%Jacob(x)%*%s+t(Jacob(x)-Jacob(x_pre))%*%f(x)

L=Jacob(x)+W;B=t(L)%*%L

W=W+(L%*%s/as.numeric(c(s)%*%B%*%c(s)))%*%t(c((sqrt(as.numeric(c(s)%*%B%*%c(s))/sum(s*z))*z-B%*%s)))

d=-solve(t(Jacob(x)+W)%*%(Jacob(x)+W))%*%t(Jacob(x))%*%f(x)

x_pre=x

x=x+eta*d

print(sum(abs(f(x))))

}


```

```{r}

#p.102~p.104 Huschens method (学習率については今後学習予定(Armijoの基準))

f=function(s){

x1=s[1];x2=s[2]

f1=x1^2-x2

f2=x1-x2^2

return(c(f1,f2))

}

Jacob=function(s){

vec=s

J=array(0,dim=c(length(s),length(s)))

for(j in 1:length(s)){

vec_sub=vec;vec_sub[j]=vec_sub[j]+h

J[,j]=(f(vec_sub)-f(vec))/h

}

return(J)

}


ite=10^2

x=rep(5,2)

x_pre=rep(0,2)

eta=0.1

W=diag(rep(1,length(x)))

h=0.01

A=diag(rep(0,length(x)))

for(l in 1:ite){

d=-solve(t(Jacob(x))%*%Jacob(x)+sum(f(x)^2)*A)%*%t(Jacob(x))%*%f(x)

x_pre=x

x=x+eta*d

s=x-x_pre

q=t(Jacob(x)-Jacob(x_pre))%*%(f(x)/sum(f(x)^2))

z=t(Jacob(x))%*%Jacob(x)%*%s+sum(f(x)^2)*q

B=t(Jacob(x))%*%Jacob(x)+sum(f(x)^2)*A

A=A+(1/sum(f(x)^2))*(t(t(c(z)))%*%t(c(z))/sum(s*z)-B%*%t(t(c(s)))%*%t(c(s))%*%B/as.numeric(c(s)%*%B%*%c(s)))

print(sum(abs(f(x))))

}


```