この記事はレイトレアドベントカレンダー2021の11日目の記事です。

逆関数無しの逆関数法

この投稿では数値積分法のモンテカルロ法の逆関数法の新しい手法[1]について書き記します。
普通の逆関数法は、原始関数に逆関数がなく直接的に逆関数法が使えないような場合使えません。
しかしいくつかの条件を満たすことで、そんな関数に対しても逆関数法を行えるようになります。

モンテカルロ法

モンテカルロレイトレーシングでは、モンテカルロ法を用いてレンダリング方程式を数値計算的に解きます。モンテカルロ法では、対象の関数の定義域を一様にサンプルすることで近似的にその関数の積分値を求めます。

$$ \bar{\theta} = \frac{1}{n} \sum^{n}_{i=1}{f(x_i)} $$

例えば次の関数を[0,1]の範囲で積分すると1になります $$ f(x) = \frac{\pi}{2 \sqrt{2}} \sin(\frac{\pi}{2} (x+\frac{1}{2})) $$

この計算をモンテカルロ法で行うと次のようになります。

In [40]:
import math
import random
import matplotlib.pyplot as plt

# 便利な定義
sin  = math.sin
cos  = math.cos
asin = math.asin
pow  = math.pow
exp  = math.exp
log  = math.log
sqrt = math.sqrt
pi   = math.pi

# 積分したい関数。[0,1]で積分すると1になる
def f_sin(x): return sin((x+0.5) * 0.5 * pi) * 0.5 * pi / sqrt(2.0)

# プロット
fig = plt.figure(figsize=(6.0, 6.0))
ax  = fig.add_subplot(1, 1, 1)
ax.set_xlim([0.0,1.0])
ax.set_ylim([0.0,1.2])
x = [x/100.0 for x in range(100)]
y = [f_sin(x) for x in x]
ax.plot(x,y, label="f(x)")

# 一様にサンプル
est = 0.0
random.seed(0x123)
sx = [random.random() for i in range(32)]
sy = [f_sin(x) for x in sx]
est = sum(sy)/len(sy)
ax.vlines(
    x=sx, 
    ymin=[0 for _ in range(len(sx))],
    ymax=[sy],
    color='red', alpha=0.5, lw=4, 
    label='estimation:{:.2f}'.format(est))
ax.legend(loc = 'upper right')
plt.show()

関数$f(x)$のプロットの下に、サンプルした位置$x_i$に、その位置での関数$f(x_i)$の値を長さとして赤い線を引きました。凡例のところに推定値が掛かれています。1なので正しい答えを得られていそうです！

しかし単に一様にサンプルするだけだと効率が悪いケースがあります。具体的にはほとんどの領域でほぼ0となるような関数では一様にサンプルすると、ほとんどのサンプルでほぼ0となってしまい効率が悪いです。例えば次のようなテント関数を考えてみます。

$$ f(x) = \begin{cases} 4x &\text{if } x<0.5 \\ 4x(1-x) &\text{if } x\geqq0.5 \end{cases} \\ $$

コードにすると次のようになります。

In [41]:
# テント関数。[0,1]で積分すると1になる
def f_tent(x): 
    if x < 0.5 : return 4.0 * x
    else:        return 4.0 * (1.0 - x)

この関数を[0,1]の範囲で積分すると1になります。この関数のプロットとモンテカルロ積分をしてみます。

In [42]:
import math
import numpy as np
import matplotlib.pyplot as plt

# プロット
fig = plt.figure(figsize=(6.0, 6.0))
ax  = fig.add_subplot(1, 1, 1)
ax.set_xlim([0.0,1.0])
ax.set_ylim([0.0,2.2])
x = np.linspace(0, 1, 100)
y = [f_tent(x) for x in x]
ax.plot(x,y, label="f(x)")

# 一様乱数でサンプル
est = 0.0
random.seed(0x1234)
sx = [random.random() for i in range(32)]
sy = [f_tent(x) for x in sx]
est = sum(sy)/len(sy)
ax.vlines(x=sx, 
          ymin=[0 for _ in range(len(sx))], 
          ymax=[sy],
          color='red', alpha=0.5, lw=4, 
          label='estimation:{:.3f}'.format(est))
ax.legend(loc = 'upper right')
plt.show()

この計算をナイーブなモンテカルロ法で行うと、値がほとんどないところにもサンプルが生成され、結果推定値が1から大きくずれるのがわかります。

重点サンプリング

このような場合、もし$f(x)$に比例する確率でサンプルできれば、一様にサンプルした場合に比べて効率がよさそうです。つまり、より$f(x)$が大きいところはより多くサンプルし、$f(x)$が小さいところはより少なくサンプルを少なくしたいです。そのサンプルする確率を$p(x)$とします。サンプルする確率を一様から変えてしまうと、$p(x)$が大きいところの結果が平均値に影響を与えすぎるので、$p(x)$で割ることで、相殺します。これを重点サンプリングと呼びます。

$$ \bar{\theta} = \frac{1}{n} \sum^{n}_{i=1}\frac{f(x_i)}{p(x_i)} $$

「関数の値に比例する確率でサンプルする」と書きましたが、それはどうやるのでしょうか？ それには棄却法や逆関数法が使えますが、ここでは逆関数法について書きます。

逆関数法

逆関数法では次のようにサンプルxを決めます。

  1. まず$f(x)$の定義域全体で積分結果が1になるような関数$p(x)$を用意します。これを確率密度関数(PDF,Probability density function)と呼びます。PDFは$f(x)$に近い形状が良いですが、必ずしも$f(x)$と同じである必要はありません。($f(x)$が0ではないところは$p(x)$も必ず0以外になっている必要もあります)
  2. このPDFを積分し原始関数を求めます。PDFの原始関数の事を累積分布関数(CDF, Cumulative distribution function)と呼びます。
  3. CDFの逆関数を求めます。
  4. CDFの逆関数に[0,1)の乱数$u$を入れることで、$p(x)$の値に比例した確率でサンプル$x$が得られます。
  5. そのままだと$p(x)$分だけサンプルし過ぎになるので、$p(x)$で割って相殺します。

文字だけだと、とてもわかりづらいので具体例で示します。

先のテント関数を再掲します。

$$ f(x) = \begin{cases} 4x &\text{if } x<0.5 \\ 4x(1-x) &\text{if } x\geqq0.5 \end{cases} \\ $$

この関数の原始関数$F(x)$とさらにその逆関数$F^{-1}(u)$は次の通りです。

$$ F(x) = \begin{cases} 2x^2 &\text{if } x<0.5 \\ 2x(2-x) &\text{if } x\geqq0.5 \end{cases} \\ F^{-1}(u) = \begin{cases} \sqrt{\frac{u}{2}} &\text{if } u<0.5 \\ 1-\sqrt{\frac{1-u}{2}} &\text{if } u\geqq0.5 \end{cases} \\ $$

コードにすると次の通りです。

In [43]:
# p(x)の原始関数の逆関数
def sample_tent(u) : 
    if u < 0.5 : return sqrt(u/2.0)
    else :       return 1.0 - sqrt((1.0-u)/2.0)

# 注意！ 今回はp(x)==f(x)の特別な例です
def p_tent(x):
    return max(f_tent(x),1e-10)

これを使ってサンプルしてプロットしてみます。

In [44]:
import math
import numpy as np
import matplotlib.pyplot as plt

# プロット
fig = plt.figure(figsize=(6.0, 6.0))
ax  = fig.add_subplot(1, 1, 1)
ax.set_xlim([0.0,1.0])
ax.set_ylim([0.0,2.2])
x = np.linspace(0, 1, 300)
y = [f_tent(x) for x in x]
ax.plot(x,y, label="f(x)")

# 重点サンプリング
est = 0.0
random.seed(0x12345)
sx = [sample_tent(random.random()) for i in range(32)]
sy = [f_tent(x) for x in sx]
est = sum([f_tent(x)/p_tent(x) for x in sx])/len(sx)
ax.vlines(x=sx, 
          ymin=[0 for _ in range(len(sx))],
          ymax=[sy], color='red', alpha=0.5, lw=4, 
          label='estimation:{:.2f}'.format(est))
ax.legend(loc = 'upper right')
plt.show()

関数の値の高いところにサンプルが集中し、なおかつ推定値が1に近づいていることが分かります。

さらに、サンプルした数が本当にこの関数の形に比例するかをヒストグラムで確認します。

In [45]:
# vdc
def ri(n, d):
    base = [2,3,5,7,11,13,17,19,23,29][d]
    nume = 0
    deno = 1
    while n > 0:
        nume = nume * base + (n % base)
        deno = deno * base
        n //= base
    return float(nume) / float(deno)
    
# 元のpdfと、それに対するサンプルのヒストグラムのプロット
def plotHist(f,samples):
    # 関数
    fig = plt.figure(figsize=(6.0, 6.0))
    ax  = fig.add_subplot(1, 1, 1)
    x = [x/100.0 for x in range(100)]
    y = [f(x) for x in x]
    ax.plot(x,y, label="f(x)")
    # ヒストグラム
    plt.hist(samples,density=True, bins=50,color='cornflowerblue', edgecolor='black')
    plt.show()

# サンプル
samples = [sample_tent(ri(i,0)) for i in range(2048)]
# ヒストグラムを表示
plotHist(f=f_tent, samples=samples)

ばっちり一致しています！

このように逆関数法はとても便利です。どのような関数でもサンプルをするには逆関数法を使えばよいのではないか？という気分になりますが、しかしここで大きな問題があります。原始関数の逆関数が求まらない関数が非常に多いのです。

欠けた半円問題

例えば、欠けた半円 のような関数について考えます。

$$ f(x)=\sqrt{1-(ax)^2} $$

この関数の原始関数は次の通りです。

$$ F(x)=\frac{x}{2} \sqrt{1-(ax)^2} + \frac{1}{2a} \sin^{-1}(ax) $$

しかしこの関数$F(x)$には逆関数は存在しません。形状としてはシンプルに感じる関数ですら意外なことに原始関数の逆関数は求まらないのです。欠けた半円の$f(x)$と$F(x)$をコードにすると次のようになります。後々使うので導関数$f'(x)$も実装しておきます。

In [46]:
import math
# 
a = 0.90
# 欠けた半円の正規化項
def F_1_td(): return (0.5 * sqrt(1.0 - a * a) + (0.5 * asin(a))/a)
# 欠けた半円(truncated disk)。[0,1]で積分すると1になる
def f_td(x): return sqrt(1.0-x*x*a*a) / F_1_td()
# f(x)の原始関数
def F_td(x): return (0.5 * x * sqrt(1.0 - a * a * x * x) + (0.5 * asin(a*x))/a) / F_1_td()
# p(x)の導関数
def fp_td(x): return -(a*a*x)/(sqrt(1.0-a*a*x*x)) / F_1_td()

この関数をプロットしてみます

In [47]:
import math
import matplotlib.pyplot as plt

COLOR_REGION_A = 'cornflowerblue'
COLOR_REGION_B = 'lightgreen'
COLOR_LINE     = 'black'

#
def plotTruncatedDisk(fillAB=False, xa=None, xb=None, labela=None, labelb=None, scatterSamples=None):
    # プロット
    width = 5.0
    fig = plt.figure(figsize=(width, width*1.2))
    ax  = fig.add_subplot(1, 1, 1)
    ax.set_xlim([0,1.0])
    ax.set_ylim([0,1.2])
    numGrid = 300
    xs = np.linspace(0, 1, numGrid)
    ys = np.array([f_td(x) for x in xs])
    ax.plot(xs,ys, label="f(x)", color=COLOR_LINE)

    # 分割する領域を示す
    if xa != None and xb != None:
        yb = 0.0
        ya = f_td(xa)
        plt.plot([xa, xb], [ya, yb], color=COLOR_LINE)
        if xa != xb:
            slope = (yb-ya)/(xb - xa)
            ssy = np.array([slope * (x-xb) for x in xs])
        else:
            ssy = np.array([(x>=xb) for x in xs])
        y2 = ys
        y1 = np.maximum(np.zeros(numGrid), ssy)
        y3 = np.zeros(numGrid)
        y4 = np.minimum(ys,ssy)
        ax.fill_between(xs, y1, y2, facecolor=COLOR_REGION_A, where= y2>=y1)
        ax.fill_between(xs, y3, y4, facecolor=COLOR_REGION_B, where= y3<=y4)

        #
        if labela != None:
            ax.text(0.3, 0.4, labela, fontsize=25)
        if labelb != None:
            ax.text(0.7, 0.4, labelb, fontsize=25)
    
    # 散布図
    if scatterSamples != None:
        X = scatterSamples[0]
        Y = scatterSamples[1]
        ax.scatter(X,Y,s=5, color=COLOR_LINE)
    # 
    plt.show()

plotTruncatedDisk()

この欠けた半円$f(x)$の原始関数の逆関数は求まらないので、直接逆関数法を使うことはできません。

逆関数法が使えないからサンプルはできないというわけではありません。原始関数$F(x)$は求まるので、これを対象にして$u=F(x)$を満たすサンプル$x$を二分探索をすることはできます。実装は次の通りです。

In [48]:
# 二分探索で F^{-1}(u) = x を求める
def findX(u):
    ok = 0.0
    ng = 1.0
    while ng - ok > 1E-10:
        mid = (ok + ng)/2.0
        if F_td(mid) < u: ok = mid
        else:             ng = mid
    return ok

u = 0.7
x = findX(u)
plotTruncatedDisk(xa=x,xb=x, labela="u", labelb="1-u")

二分探索することで、$f(x)$が$u:1-u$の領域に確かに分割されました。このときこの区切り線のx座標がそのままサンプルになります。しかし、二分探索は計算コストが高いので解析的にサンプルできるとよりうれしいです。Triangle-Cut Parameterizationはこの逆関数法が使えないような関数に対しても別の角度から問題を解決し、解析的にサンプルできるようにします。

二つの領域を分ける

まず$f(x)$を直接使うのではなく、逆関数法を使えるよりシンプルな関数$g(x)$を導入します。 これは$f(x)$により近い形状がうれしいのですが、とりあえずここでは$g(x)=1$とします。この$g(x)$を使って逆関数法を使うために、その原始関数$G(x)$、その逆関数$G^{-1}(u)$を用意します。

In [49]:
# 逆関数法が使えるシンプルな関数
def g_const(x) : return 1.0
# g(x)の原始関数
def G_const(x) : return x
# G(x)の逆関数
def iG_const(u) : return u

$f(x)$の事はいったん脇に置いて、$g(x)$を使ってサンプルしてしまいます。すると次のようなプロットになります。

In [50]:
# 近似のXを求めてプロット
u = 0.7
xa = iG_const(u)
plotTruncatedDisk(xa=xa,xb=xa, labela="iG(u)", labelb="1-iG(u)")

当然、この左右の領域の面積の比は$u:1-u$ではありません。具体的にずれた誤差は、 $$ err = u - F(x_a) \\ where \: x_a = G^{-1}(u) $$ となります。ここで、区切り線を斜めにすることを考えます。$x_a$と$F(x)$の交点は固定して、$x_a$とX軸との交点、つまり切片の方を$x_b$として動かすことを考えます。その時にできる三角形の面積は次のように表されます。 $$ area = \frac{1}{2} (x_b - x_a) f(x_a) $$ この三角形の面積と誤差が一致しているとき、領域の比は$u:1-u$となります。そのときの$x_b$は次のようになります。 $$ x_b = x_a + 2\frac{u - F(x_a)}{f(x_a)} $$ とても簡単な式になりました。これを使ってプロットしてみます。

In [51]:
# X切片を求める
u = 0.7
xb = xa + 2.0 * (u-F_td(xa)) / f_td(xa)
# プロット
plotTruncatedDisk(xa=xa,xb=xb, labela="u", labelb="1-u")

この線によって、領域が$u:1-u$に確かに分割されました！

しかし問題があります。いつもの逆関数法であればxの値が求まって終わりなのですが、今回は区切りの線が斜めです。この線のどこをサンプルすればいいのかということがわかりません。ここからさらにこの区切りの線の上をサンプルすることを考えます。

無限小の線分の上のサンプル

分割線の上をサンプルします。この分割線を$u$を微小量$du$だけ増やした時に、点$a$と点$b$が動いた量を$w_a$、$w_b$と置きます。この時掃いた部分には、無限小の幅を持つ台形ができます。この時の$w_a$と$w_b$の関係は天下りですが次のようになります。

$$ \begin{align} & w_a \propto f(x_a)^2 + 2(u-F(x_a)) f'(x_a) \\ & w_b \propto 2f(x_a)g(x_a) - w_a \tag{1} \end{align} $$

この無限小の幅を持つ台形をサンプリングしていきます。 $w_a$と$w_b$の比例関係をそのまま等式だと思うと、その幅の大きさ$w(t)$は次のようになります。

$$ w(t) = (w_b-w_a)t + w_a $$

この関数をpdfにするために、$w_a+w_b$で割ってしまいます。

$$ w(t) = \frac{w_b-w_a}{w_a+w_b}t + \frac{w_a}{w_a+w_b} $$

これで$w(t)$はpdfとなりました。さらにこれの原始関数は次のようになります。

$$ W(t) = \frac{w_b-w_a}{2(w_a+w_b)}t^2 + \frac{w_a}{w_a+w_b} t $$

$F(x)$の逆関数は二次関数の解の的な形になりますが、数値計算的に安定させるために式変形の途中で有理化させます。結果は次の通りになります。

$$ W^{-1}(v) = \frac{v(w_a+w_b)}{w_b + \sqrt{(1-v)w_b^2 + vw_a^2}} $$

$w(t)$に関する関数一式が揃いました。百聞は一見に如かず、この関数を使ってサンプルしてみましょう！

In [52]:
# 台形の端点の厚さ
wa = 0.1
wb = 0.5
# 台形の厚さの関数
def w_trapezoid(x): return (wb-wa)*x + wa
# w(x)の原始関数の逆関数
def iW_trapezoid(u): return u*(wa+wb)/(wa+sqrt((1-u)*wa*wa + u*wb*wb))

# プロット
fig = plt.figure(figsize=(6.0, 6.0))
ax  = fig.add_subplot(1, 1, 1)
ax.set_xlim([0.0,1.0])
ax.set_ylim([0.0,0.5])
x = np.linspace(0.0, 1.0, 100)
y = [w_trapezoid(x) for x in x]
ax.plot(x,y, label="f(x)")

# サンプル
sx = [iW_trapezoid(ri(i,0)) for i in range(32)]
sy = [w_trapezoid(x) for x in sx]
est = 0.0
ax.vlines(x=sx, 
          ymin=[0 for _ in range(len(sx))],
          ymax=[sy], color='red', alpha=0.5, lw=4, 
          label='estimation:{:.2f}'.format(est))

plt.show()

台形の高さに応じてサンプルできている様子が確認できました。
これで線分上を適切にサンプルできるようになりました。

これによって、$f(x)$の中を一様にサンプルすることができるようになりましたが、最初欲しかったものは$F^{-1}(u)$でありスカラー値でした。どうやって平面上のサンプルを一次元にすればいいのでしょうか。実はこれは簡単で、単にX座標を見るだけで大丈夫です。

長かったですが、これですべてのパーツが出そろいました！
いよいよ一つの仕組みに落としていきます。

Triangle-Cut Parameterization

まとめると次のような話の流れでした。

  1. $f(x)$に従う確率でサンプルしたい
  2. 二つの乱数$(u,v)$を用意する
  3. 原始関数とその逆関数が求まるシンプルな関数$g(x)$を用意する。
  4. $G(x)$と$F(x)$を使って$f(x)$を$u:1-u$に2分割する
  5. 分割線の上を微小な面の台形とみなして、$v$を使ってサンプルする。この時$f'(x)$を用いる。
  6. $(u,v)$から$f(x)$の下に一様に分布する2次元上の座標が得られたので、その座標のX軸をサンプルとする。

これまで処理を一つの関数にまとめます

In [53]:
import math
# 
def trianglecut(u, v, f, F, fp, g, iG):
    # 点Aの座標を計算する
    xa     = iG(u)
    f_xa   = f(xa)
    F_xa   = F(xa)
    fp_xa  = fp(xa)
    ya     = f_xa
    g_x_a  = g(xa)
    # 分割する線の座標を求める
    xb = xa + 2.0 * (u - F_xa) / f_xa
    yb = 0.0
    # 微小な分割線の点Aと点Bでの厚みを出す
    wa = f_xa * f_xa + 2.0 * (u- F_xa) * fp_xa
    wb = 2.0 * f_xa * g_x_a - wa
    # vに沿ってサンプル
    #print(wa,wb,v)
    t = v * (wa+wb) /(wb + sqrt((1.0-v) * wb * wb + v * wa * wa))
    # 最終的な座標を出す
    x = xa * t + xb * (1.0 - t)
    y = ya * t + yb * (1.0 - t)
    #
    return [x,y]

一つの関数にまとめるとかなりシンプルになりました！
ちゃんと一様に分布しているかを確認するために散布図をプロットしてみます。

In [54]:
#  Triangle-Cut Parameterizationを使ったサンプル
def sampleWithTriangleCut(numSample, f, F, fp, g, iG):
    X = []
    Y = []
    for si in range(numSample):
        u = ri(si,0)
        v = ri(si,1)
        [x,y] = trianglecut(u=u,v=v, f=f, F=F, fp=fp, g=g, iG=iG)
        X.append(x)
        Y.append(y)
    return [X,Y]

# サンプルをプロット
samples = sampleWithTriangleCut(1024, f=f_td, F=F_td, fp=fp_td, g=g_const, iG=iG_const)
plotTruncatedDisk(scatterSamples=samples)

無事一様に分布していることを確かめられました！
最後に、当初やりたかったX軸に沿ってのサンプルした結果のヒストグラムを見てみましょう。

In [55]:
# サンプル
samples = sampleWithTriangleCut(8192, f=f_td, F=F_td, fp=fp_td, g=g_const, iG=iG_const)
# ヒストグラムを表示
plotHist(f=f_td, samples=samples[0])

ピッタリと$f(x)$に比例したサンプルできていることがわかりました！

念のため別の関数でも動作を確かめてみます。

In [56]:
# 二次関数
def f_quad(x): return (x*x + x + 1.0)*6.0/11.0
def F_quad(x): return (x*x*x/3.0 + x*x/2.0 + x)*6.0/11.0
def fp_quad(x): return (2.0*x + 1)*6.0/11.0
# サンプル
samples = sampleWithTriangleCut(8192, f=f_quad, F=F_quad, fp=fp_quad, g=g_const, iG=iG_const)
# ヒストグラムを表示
plotHist(f=f_quad, samples=samples[0])

大丈夫そうです。

制限事項

無事、円が欠けたケースではサンプルできたのですが、実は$g(x)$の選び方がうまくないと、この手法は正しく動作しません。

ダメなケース1

$(u,v)$がTriangle-Cutによって$(x,y)$にマップされるわけですが、この領域が必ずしも$f(x)$の定義域に入る保証はありません。全ての$(u,v)$が出す$(x,y)$対して、$y \leq f(x)$が成り立っている事を確認する必要があります。

(機械的にはできずに、手あたり次第$(u,v)$をマップする以外に確認方法はないのだろうか...?)

ダメなケース2

無限小の幅を持つ台形について述べていたところの議論は、全て$w_a$,$w_b$が常に正であることを前提にしていました。つまり$u$の増加に伴い常に右側に移動するようなものでないといけませんでした。つまり$(1)$式が常に正である必要があります。もしこの前提が崩れるとサンプルされない領域が発生します。

二次元以上への拡張

これまでの話は1次元の関数に対してだけでしたが、2次元以上へも拡張が可能です。詳細は論文[1]を参照してください。

おわり

長い文章になってしまいましたが、コアのアイデアはシンプルで、使われる算数も簡単なモノだけで完結しているのに、逆関数法を拡張できたのは純粋に驚きでした。いくつかの関数を実際に試した肌感覚では、制限事項は見た目よりもシビアな印象でしたが、せっかく新しい道具が手に入ったのですから、ぜひいろんな関数のサンプルをしてみましょう！

リファレンス

[1] Heitz, E. (2020). Can’t Invert the CDF? The Triangle-Cut Parameterization of the Region under the Curve. Computer Graphics Forum, 39(4), 121–132. https://doi.org/10.1111/cgf.14058