optim関数で数値的に逆関数の計算したい

これがオーソドックスな正しいやり方なのかわからないけど、とりあえずやってみた。

optim関数で数値的に逆関数の計算したい(一次元版)

とりあえず一次元の関数 $f:\mathbb{R} \rightarrow \mathbb{R}$ を $y = f(x)$ と書いた時の、 $y$ が与えられた下での $x$ の推定を考えたい。

これはもちろん解析解があるとハッピーなわけだが、ない場合には数値的にやるしかないわけです。ここでは適当に $f(x) = \sin(x) + 2\sin(x)^2$ って関数を題材にやってみんとす。
この関数 $f$ をRで書くと以下のような感じか。

#適当な関数f
f <- function(x)
{
  A <- 1:2
  A[1]*sin(x) + A[2]*sin(x)^2
}

描画すると、こんな感じ。

a <- seq(-5, 5, 0.1)
plot(a, sapply(a, f))

とりあえずこの関数 $f$ に対する逆関数を計算してみたい。

そこで
$\hat{x} = arg\,min_{x} \left( f(x) - y \right)^2$
となるような $\hat{x}$ を計算するコードをoptim関数*1を使って以下のように書いた。要するに「最小二乗法でそれっぽい $x$ 探してみようや！」ですわ。
とりあえず真の $x$ は2と設定し、 $y=f(x=2)$ を求めた上で、その $y$ を使って $x$ を計算する。

#目的関数（二乗誤差）の生成
make.objective <- function(y)
{
  function(x){
    (f(x)-y)^2
  }
}
#真のx=2.0
x.answer <- 2.0
y <- f(x.answer)
#与えられたyに対してxを推定（最小二乗誤差を計算、xの初期値1.6と設定）
opt <- optim(c(1.6), make.objective(y), method="BFGS")

この最適化結果を見てみると

> #推定したパラメータ値x
> opt$par
[1] 1.999999
> #最適化した結果（二乗誤差）
> opt$value
[1] 2.932982e-12

となっていて、 $x$ の推定値であるopt$parが1.999となっているので、真の値2に近いですね、合ってそうですねと理解できる。

一方、初期条件のみを変えてもう一回実行してみると…

> #与えられたyに対してxを推定（最小二乗誤差を計算、xの初期値2.5と設定）
> opt <- optim(c(2.5), make.objective(y), method="BFGS")
> #推定したパラメータ値x
> opt$par
[1] -4.283186
> #最適化した結果（二乗誤差）
> opt$value
[1] 3.104426e-12

$x$ の推定値であるopt$parが-4.28程度になってしまった！これは上のグラフを見てもわかるように $y=2$ となる点が $x=-4.28$ 辺りも含めて複数あるからで、まぁ、これはどうしようもない気もするが、世の中的にはどうやるのが普通なのか*2、よく知らねぇわ…。

optim関数で逆関数の計算したい(多次元版)

先ほどの関数 $f$ を二次元に拡張して同じようなことをやって見る。
目的関数を多次元の関数 $f:\mathbb{R}^2 \rightarrow \mathbb{R}^2$ として、さきほど同様、 $y = f(x)$ と書いた時の、 $y$ が与えられた下での $x$ の推定を考えたい。関数fはさっきの亜種で以下のように定義。

f <- function(x)
{
  A <- 1:4
  y1 <- A[1]*sin(x[1]) + A[2]*sin(x[1])^2
  y2 <- A[3]*sin(x[2]) + A[4]*sin(x[2])^2
  c(y1, y2)
}

目的関数、正確にはその生成関数も、全体での二乗誤差が最小になるように変更しておく。

#目的関数（二乗誤差）の生成
make.objective <- function(y)
{
  function(x){
    sum((f(x)-y)^2)
  }
}

とりあえず真の $x=(2, 2.5)$ と設定し、[tex:y=f*3]を求めた上で、その $y$ を使って $x$ を計算する。

#真のx=c(2.0, 2.5)
x.answer <- c(2.0, 2.5)
y <- f(x.answer)
#与えられたyに対してxを推定（最小二乗誤差を計算）
opt <- optim(c(1.8, 2.4), make.objective(y), method="BFGS")

で結果は

> opt$par
[1] 1.999999 2.500000
> #最適化した結果（二乗誤差）
> opt$value
[1] 2.948609e-12

たしかにopt$parが $x=(2, 2.5)$ に近いので良さそうですよと、まぁ、これはよしと。で、初期値をちょろっとずらすと

> opt <- optim(c(1, 0.3), make.objective(y), method="BFGS")
> #推定したパラメータ値x
> opt$par
[1] 1.9999991 0.6415926
> #最適化した結果（二乗誤差）
> opt$value
[1] 3.160292e-12

あぁ、先ほど同様、にだめだねと。理由も一次元の場合同様だろう。

念の為、最小二乗誤差を可視化する。
最小二乗誤差は以下のようにggplot2を使うと等高線＋濃淡(ヒートマップ？)で可視化できて、真の $x$ の箇所以外にも、ダメったケース(2,0.64)も含めて、誤差が小さくでる箇所（背景色濃）の存在を確認できる。

a <- seq(-3, 3, 0.1)
obj <- make.objective(y)
grid <- expand.grid(x=a, y=a)
df <- cbind(grid, apply(grid, 1, obj))
names(df) <- c("x", "y", "z")
library(ggplot2)
ggplot(df, aes(x,y,z=z))  + 
  geom_tile(aes(fill=z)) +
  stat_contour()+ 
  xlim(-3, 3) + 
  ylim(-3, 3)