射影定理と関数近似

記事の内容


関数解析の重要な結果の1つである射影定理と, その応用としての近似理論について簡単にまとめておきます. 最初の節で, 射影定理を導入し, 続く2節でその応用を数値実験結果とともに紹介します.

射影定理

Banach空間とHilbert空間

関数解析学は, Banach空間とHilbert空間を舞台に展開されます. 本記事では主にHilbert空間を扱います. まずはこれらの空間を定義しておきましょう.

【定義 (Banach空間とHilbert空間)】

以下の(1)から(3)を定める.

      (1) 実数体上のノルム空間\((V,\|\cdot\|)\)の任意のCauchy列が収束するとき, \((V,\|\cdot\|)\)は完備であるという.
      (2) 完備なノルム空間をBanach空間という.
      (3) 実数体上の内積空間\((\mathcal{H},(\cdot,\cdot))\)は, 内積によって誘導されるノルム\(\|\cdot\|_{\mathcal{H}}=\sqrt{(\cdot,\cdot)}\)に関して完備であるとき, Hilbert空間であるという.

以下, \(L^2([a,b])\)上の内積を以下のように定めます:

\begin{equation} (f,g)_2 = \int_{[a,b]} f(x) g(x) \mu(dx). \end{equation} ただし, \(\mu\)はLebesgue測度です. このとき, \((L^2([a,b]),(\cdot,\cdot)_2)\)はHilbert空間です.

射影定理

次に, Hilbert空間における射影定理を述べておきます.

【定理 (射影定理)】

Hilbert空間\((\mathcal{H},(\cdot,\cdot))\)の閉部分空間\(M\)について, 任意の\(y\in\mathcal{H}\)に対して, ある\(\hat{y}\)が存在して, 次式が成り立つ: \begin{equation} \hat{y} = \underset{m\in M}{\mathrm{argmin}}{\|y-m\|}. \end{equation} さらに, 任意の\(m\in M\)に対して\((y-\hat{y},m)=0\)が成り立つことと, \(y\)に対する\(\hat{y}\)の存在の一意性は同値である.

部分空間の閉性は以下のように保証されます.

【命題 (有限次元部分空間は閉である)】

ベクトル空間の有限次元部分空間は閉である.

応用1:最小二乗近似

この節では, 射影定理の応用として, 最小二乗近似を説明します. 機械学習系の書籍では正規方程式を微分により導出しますが, ここでは射影定理との関係を強調した形で述べます.

問題設定と正規方程式の導出

問題設定を確認しておきましょう. Hilbert空間\((\mathcal{H},(\cdot,\cdot)\)の内積から誘導されるノルムを\(\|\cdot\|_{\mathcal{H}}\)と表しておきます. この空間の既知の元\(\phi_1,\dots,\phi_r\)によって生成される部分空間を\(M\)と表します. このとき, 任意の\(y\in\mathcal{H}\)を既知のベクトル\(\phi_1,\dots,\phi_r\)の一次結合で近似する問題, つまり, 次の最小化問題を考えます:

\begin{equation} w_1,\dots,w_r = \underset{w_1,\dots,w_r}{\mathrm{argmin}} \| y-w_1\phi_1-\cdots-w_r\phi_r\|_{\mathcal{H}}. \end{equation}

部分空間\(M\)は有限次元なので, 閉であり, 射影定理により上の最適化問題の解の存在が保証されます. さらに, 解の存在が一意ならば, 各\(j=1,\dots,r\)に対して次式が成り立ちます:

\begin{equation} 0 = (y-w_1\phi_1-\cdots-w_r\phi_r,\phi_j) = (y,\phi_j) - w_1(\phi_1,\phi_j)-\cdots-w_r(\phi_r,\phi_j). \end{equation}

したがって, 次式が成り立ちます:

\begin{equation} \begin{bmatrix} (\phi_1,\phi_1) & \cdots & (\phi_r,\phi_1) \\ \vdots & \ddots & \vdots \\ (\phi_1,\phi_r) & \cdots & (\phi_r,\phi_r) \end{bmatrix} \begin{bmatrix} w_1\\ \vdots\\ w_r \end{bmatrix} = \begin{bmatrix} (y,\phi_1)\\ \vdots\\ (y,\phi_r) \end{bmatrix}. \end{equation}

この方程式は正規方程式と呼ばれ, この解が求めたい係数です. 空間とノルムを適切に定義することで, 様々な問題に適用でき, ここでは, いわゆる"最小二乗法"と, 最小二乗多項式の問題を考えます. ちなみに, 基底ベクトル\(\phi_1,\dots,\phi_r\)の一次独立性と, 解の一意性は同値です.

例1:いわゆる最小二乗法

いわゆる"最小二乗法"を考えます. 既知の入力データ\(x_1,\dots,x_N\in\mathbb{R}^{D}\)に対して, 既知の基底関数\(\phi_1,\dots,\phi_r\)を通して, 既知の出力データ\(y=[y_1,\dots,y_N]^{\mathrm{T}}\in\mathbb{R}^N\)によく当てはまる直線を計算します. すなわち, 標準内積の入ったHilbert空間\((\mathbb{R}^N,(\cdot,\cdot))\)において, 二乗距離の意味で以下のように近似します:

\begin{equation} \begin{bmatrix} y_1\\ \vdots \\ y_N \end{bmatrix} \simeq w_1 \begin{bmatrix} \phi_1(x_1)\\ \vdots \\ \phi_1(x_N) \end{bmatrix} + \cdots + w_r \begin{bmatrix} \phi_r(x_1)\\ \vdots \\ \phi_r(x_N) \end{bmatrix} = \Phi w. \end{equation}

ただし, \(\Phi\)は以下のように定めます:

\begin{equation} \Phi = \begin{bmatrix} \phi_1(x_1) & \cdots & \phi_r(x_1) \\ \vdots & \ddots & \vdots \\ \phi_1(x_N) & \cdots & \phi_r(x_N) \end{bmatrix}. \end{equation}

これは計画行列と呼ばれ, 正規方程式は以下のように表示できます:

\begin{equation} \Phi^{\mathrm{T}}\Phi w = \Phi^{\mathrm{T}}y. \end{equation}

実験してみます. サイズが\(N=25\)の入力データ\(\{x_i\}_{i=1}^{N}\)を-3から3の間の一様乱数, 出力データを次の式から生成します:

\begin{equation} 0.7\{ \sin(3x) + \cos(2x) + \sin(5x)\} + \varepsilon,\quad \varepsilon\sim\mathrm{N}(0,0.5^2). \end{equation}

上式は本来未知であり, 近似の対象です. 基底関数を以下のように定めます:

\begin{equation} \phi_j(x) = \exp\left\{ -\frac{1}{2}(x-\mu_j)^2\right\},\quad \mu_j = j-4\quad (j=1,\dots,7). \end{equation}

近似の結果が以下のとおりです. 真の関数(赤実線), 最小二乗近似関数(緑破線), データ(青バツ)を示しました. 基底関数の配置にも依ますが, 結構外していますね.

【コード3の実行結果】

例2:最小二乗多項式

有界閉区間上の関数を, 有限次数の多項式で近似する問題を考えます. いま, Hilbert空間\((L^2([0,1]),(\cdot,\cdot)_2)\)上の関数列 \begin{equation} \phi_j = x^{j-1},\quad (j=1,\dots,r+1) \end{equation} をとり, 関数\(f\in L([0,1])\)を, 高々r次の多項式で, この内積が誘導するノルムに関して最小二乗近似してみましょう. すなわち, \begin{equation} f(x) \simeq w_1\phi_1 + \cdots + w_{r+1}\phi_{r+1} = w_1 + w_2x + w_3x^2 + \cdots + w_{r+1}x^{r} \end{equation} のように近似します.

正規方程式は以下の通りです:

\begin{equation} \begin{bmatrix} \frac{1}{1} & \frac{1}{2} & \cdots & \frac{1}{r+1} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{1}{r+1} & \frac{1}{r+2} & \cdots & \frac{1}{2r+1} \end{bmatrix} \begin{bmatrix} w_1\\ \vdots\\ w_{r+1} \end{bmatrix} = \begin{bmatrix} (f,\phi_1)_2\\ \vdots\\ (f,\phi_{r+1})_2 \end{bmatrix}. \end{equation}

係数行列はHilbert行列と呼ばれます. 例として, 対数関数\(\log x\)の多項式近似を考えてみましょう. 最高次数2次, 4次, 6次の場合を下図に図示しました.

【コード5の実行結果】

応用2:直交関数展開

一般化Fourier級数展開

最小二乗近似では, 正規方程式を解く必要がありますが, 一般に, 係数行列の性質は良くないので, 正規方程式は回避したいところです. 係数行列にはベクトルの内積が登場するので, 正規直交系をうまく使えば正規方程式を解かずに済みます. この方針を実現する方法として, 以下のような結果が知られています.

【定理 (一般化Fourier級数展開)】

(可分な)Hilbert空間\((\mathcal{H},(\cdot,\cdot))\)の正規直交系\(\{\phi_j\}_{j=1}^{\infty}\)に直交するベクトルがゼロベクトルのみのとき, 任意の\(y\in\mathcal{H}\)は, \begin{equation} y = \sum_{j=1}^\infty (y,\phi_j)\phi_j \end{equation} と一意に展開できる.

コンピュータで実装する際には, この和を有限で打ち切る必要があります. また, 正規直交系は, Gram-Schmidtの正規直交化法を使って構成できます.

例1:Fourier級数展開

Hilbert空間\((L^2([0,2\pi]),(\cdot,\cdot)_2)\)上の完全正規直交系 \begin{equation} \phi_j(x) = \frac{1}{\sqrt{2\pi}}\exp(\sqrt{-1}jx),\quad (j\in\mathbb{Z}) \end{equation} に関する一般化Fourier級数展開は, いわゆる"(複素)Fourier級数展開"です.

例2:Legendre多項式

Hilbert空間\((L^2([-1,1]),(\cdot,\cdot)_2)\)上の完全正規直交系 \begin{equation} \phi_j(x) = \sqrt{\frac{2j+1}{2}}\left\{ \frac{1}{2^j j!}\frac{\mathrm{d}^j}{\mathrm{d}x^j}(x^2-1)^j \right\},\quad (j=0,1,\dots) \end{equation} は, 単項式列\(\{1,x,x^2,\dots\}\)にGram-Schmidtの正規直交化を施して得られ, Legendre多項式と呼ばれます. \(\cos\)関数を0次, 2次まで近似した結果を下図に示しました. (多分)次式のように近似できます:

\begin{equation} \cos(x) \simeq \sin(1)+\frac{5}{2}(3\cos(1)-2\sin(1))(3x^2-1). \end{equation}
【コード7の実行結果】

例3:Chebyshev近似

Hilbert空間\((L^2([-1,1]),(\cdot,\cdot)_C)\)上の直交系として, Chebyshev多項式を考えます. ここで, 内積を以下のように定めます:

\begin{equation} (f,g)_C = \int f(x)g(x) \frac{\mu(dx)}{\sqrt{1-x^2}}. \end{equation}

Chebyshev多項式は以下の漸化式により定まります:

\begin{equation} T_0(x)=1,\quad T_1(x)=x,\quad T_j(x)=2xT_{j-1}(x)-T_{j-2}(x). \end{equation}

この直交系を用いて, \(f\in L^2([-1,1])\)次のように展開できます:

\begin{equation} f(x) = \frac{1}{\pi}(f,T_0)_C T_0(x) + \frac{2}{\pi}\sum_{j=1}^\infty (f,T_j)_CT_j(x). \end{equation}

関数\(x\sqrt{1-x^2}\)を1次, 3次, 5次まで展開した結果を下図に示しました. (多分)次式のように近似できます:

\begin{equation} x\sqrt{1-x^2} \simeq \frac{2}{\pi}\left\{ \frac{2}{3}x-\frac{2}{5}(4x^3-3x)-\frac{2}{21}(16x^5-20x^3+5x) \right\}. \end{equation}
【コード9の実行結果】

このChebyshev近似は, 最大値ノルム\(\|\cdot\|=\underset{x\in [a,b]}{\mathrm{max}}|\cdot|\)に関しても収束します. 最大値ノルムに関する最良近似をミニマックス近似といい, Weierstrassの多項式近似定理により, ミニマックス近似多項式の存在が保証されます. しかし, 実際にミニマックス近似多項式を求めるのは困難で, 代わりにChebyshev近似が利用されることがあります. ちなみに, 最大値ノルムはいかなる内積からも誘導されないため, Hilbert空間の枠組みでは扱うことができません.

コード

【Juliaコード1; 初期化】
#mathematics
using LinearAlgebra
using Polynomials

#statistics
using Random
using Statistics

#visualize
using Plots
pyplot()

#macros
using UnPack
using ProgressMeter
【Juliaコード2; 関数定義】
#true function
function true_func(x)
    0.7*(sin(3*x) + cos(2*x) + sin(5*x))
end

#basis function
function design_matrix(X,μs,σs,N,D)
    [exp(-0.5*(X[n]-μs[j])^2/σs[j]^2) for n in 1:N, j in 1:D]
end

#solve regular equations
function solve_regular_equations(data,params)
    @unpack X,Y,N = data
    @unpack μs,σs = params
    Φ = design_matrix(X,μs,σs,N,D)
    (Φ'*Φ)\(Φ'*Y)
end

#regression function
function reg_func(x,wvec,params)
    @unpack μs,σs = params
    dot(wvec,[exp(-0.5*(x-μs[j])^2/σs[j]^2) for j in 1:D])
end
【Juliaコード3; 実験】
#create data
Random.seed!(42)
N = 25
X = 6*(rand(N).-0.5)
Y = true_func.(X)+ 0.5*randn(N)
data = (X=X,Y=Y,N=N)

#set the parameters
D = 7
params = (μs=-3:3,σs=ones(D),D=D)

#compute approximation function 
coeff = solve_regular_equations(data,params)

#visualize the resutls
xs = -3:0.01:3
fig1 = plot(X,Y,st=:scatter,xlabel="x",ylabel="y",label="data",title="Least square",legend=:topleft,
    markershape=:x, markerstrokewidth=1, markersize=10, markercolor=:blue)
plot!(xs,true_func.(xs),label="true",color=:red,lw=2)
plot!(xs,x->reg_func(x,coeff,params),label="least square",ls=:dash,lw=2)
savefig(fig1,"figs-FA/fig1.png")
【Juliaコード4; 関数定義】
#true function
function true_func(x)
    log(x)
end

#compute Hilbert matrix
function Hilbert_matrix(max_dim)
    [1/(i+j+1) for i in 0:max_dim, j in 0:max_dim]
end

#project true function to the basis
function project_f(max_dim)
    -[1/(j+1)^2 for j in 0:max_dim]
end

#compute coefficient
function solve_regular_equations(max_dim)
    GH = Hilbert_matrix(max_dim)
    b = project_f(max_dim)
    GH\b
end

#approximate polynomial
function poly_func(x,coeffs,max_dim)
    dot([x^k for k in 0:max_dim],coeffs)
end
【Juliaコード5; 実験】
#plot the ture and approximated function
xs = 0:0.01:1
fig2 = plot(xlabel="x",title="Least square polynomial approximation",legend=:bottomright)
plot!(xs,true_func.(xs),label="true",color=:red,lw=2)
@showprogress for max_dim in 2:2:6
    coeff = solve_regular_equations(max_dim)
    plot!(xs,x->poly_func(x,coeff,max_dim),label="max_dim=$(max_dim)",ls=:dash,lw=2)
end
savefig(fig2,"figs-FA/fig2.png")
【Juliaコード6; 関数定義】
#true function 
function true_func(x)
    cos(x)
end

#Legendre polynomial expansion(0 order)
function Legendre0(x)
    sin(1)
end

#Legendre polynomial expansion(2 order)
function Legendre2(x)
    sin(1) + 5*(3*cos(1)-2*sin(1))*(3*x^2-1)/2
end
【Juliaコード7; 実験】
xs = -1:0.01:1
fig3 = plot(xs,x->true_func(x),title="Legendre approximation",xlabel="x",label="true",color=:red,lw=2)
plot!(Legendre0,label="Legendre(0 order)",ls=:dash,lw=2)
plot!(Legendre2,label="Legendre(2 order)",ls=:dash,lw=2)
savefig(fig3,"figs-FA/fig3.png")
【Juliaコード8; 関数定義】
#true function 
function true_func(x)
    sqrt(1-x^2)*x
end

#Chebyshev approximation(1 order)
function Chebyshev1(x)
    2/π*(
        2/3*x
    )
end

#Chebyshev approximation(3 order)
function Chebyshev3(x)
    2/π*(
        2/3*x + 
        (-2/5)*(4*x^3-3*x)
    )
end

#Chebyshev approximation(5 order)
function Chebyshev5(x)
    2/π*(
        2/3*x + 
        (-2/5)*(4*x^3-3*x) +
        (-2/21)*(16*x^5-20*x^3+5*x)
    )
end
【Juliaコード9; 実験】
xs = -1:0.01:1
fig4 = plot(xs,x->true_func(x),title="Chebyshev approximation",
    xlabel="x",label="true",color=:red,lw=2,legend=:bottomright)
plot!(Chebyshev1,label="Chebyshev(1 order)",ls=:dash,lw=2)
plot!(Chebyshev3,label="Chebyshev(3 order)",ls=:dash,lw=2)
plot!(Chebyshev5,label="Chebyshev(5 order)",ls=:dash,lw=2)
savefig(fig4,"figs-FA/fig4.png")
参考文献

      [1]Luenberger, 関数解析による最適理論, コロナ社, 1973
      [2]森正武, 数値解析, 共立出版, 2018
      [3]杉原正顕, 室田和夫, 数値計算法の数理, 岩波書店, 1998
      [4]Cheney, 近似理論入門, 共立出版, 1977