チェビチェフ補間 関数近似 (3)

電気

303/

303

(9)

関数近似 (3)

チェビチェフ補間

•

⊲

2012

⊲

⊲

2

•

[a, b]

f

p

f(x)

p(x)

• k · k

k f − p k

最良多項式と

Chebyshev

(2)

• P _n

n

p ^∗ = arg min _{p∈P n} k f − p k

f

k · k

n

•

arg min(arg max)

(最大)

• p ^∗ = arg min _{p∈P n} k f − p k

, k f − p k

.

,

,

として動くのは多項式

p

(

p

n

P n

).

• arg min

arg max

,

.

たとえば

http://reference.wolfram.com/language/ref/ArgMin.ja.html

最良多項式と

Chebyshev

(4)

• kk ^∞

n

n

Chebychev

• Chebychev

,

記が一定でないため

, Chebychev

.

• l 2

•

[a, b](ただし b > a)

f

n

最小二乗近似

(2)

•

{ 1, x, . . . , x ⁿ }

f

p ( x ) = P n

k =0 a _k x ^n−k

{ a 0 , a 1 , . . . , a _n }

•

•

n

P n

{ 1, x, . . . , x ⁿ }

•

f

g

(f, g) = R b

a f(x)g(x)dx

(第 7

最小二乗近似

(4)

• g kl = (x ^k , x ^l ), h k = (f, x ^k )

.

• G = (g kl ) 0 ≤k≤n,0 ≤l≤n

g kl

(n +1)

.

p

, R b

a p(x) ² dx >

0

, G

.

G

(x) =









 1

x 1

.. . . ..

x

1



















1 

















1 x . . . x

1

. ..

1











,

a = (a

, . . . , a

)

6= 0, p(x) = P

i=0

a

x

, a 6= b

R

a

a

G

(x)adt = R

a

(p(x))

dx > 0. R

a

G

(x)dt = G

, G

.

• k f − p k ^{2 2} = (f − p, f − p)

, k f − p k ^{2 2} = (f, f ) − 2 P n

k=0 a k h k + P n k=0

P n

l=0 a k a l g kl

.

• h = (h 0 , . . . , h n ) ^T , α = (a 0 , . . . , a n ) ^T

,

.

k f − p k ^{2 2} = (f, f ) − 2h ^T α + α ^T Gα.

• G

, ^∂kf−pk _∂ α ^{2 2} = 0

α

,

.

最小二乗近似

(6)

•

− h ^T + α ^T G = 0

, G

, Gα = h

.

(

).

•

,

.

•

,

,

{ 1, x, . . . , x ⁿ }

G

, a, b

,

.

•

,

.

• { r 0 (x), . . . , r _n (x) }

P _n

(

n

)

.

最小二乗近似

(8)

以下の典拠は伊理,線形代数汎論,朝倉書店, 2009および杉原,室田,線形計算 の数理,岩波書店, 2009.

•

A

k A kk A ⁻¹ k

A

•

l 2

A

A

•

•

A

A ^H A

A ^H

A

Hermite

• { r 0 (x), . . . , r n (x) }

(f, g)

, G

,

.

•

, { r 0 (x), . . . , r n (x) }

(f, g)

, ∀ k, (r k , r k ) >

0

, G

,

.

最小二乗近似

(10)

•

, (f, g)

,

(f, g) w

.

• (f, g) w = R b

a f (x)g(x)w(x)dx

.

• w(x)

.

,

, R b

a w(x)dx

,

.

• k f − p k ^w = p

(f − p, f − p) w

.

•

[a, b]

f (x)

p(x) = P n

k=0 a k r k (x)

.

, k f − p k ^2,w

p

.

•

{ r 0 (x), . . . , r n (x) }

.

最小二乗近似

(12)

• g _kl = (r _k , r _l ) _w , h _k = (f, r _k ) _w

,

G = (g kl ) 0 ≤k≤n,0 ≤l≤n

, g kl

(n + 1)

.

, h = (h 0 , . . . , h n ) ^T , α = (a 0 , . . . , a _n ) ^T

.

•

{ 1, x, . . . , x ⁿ }

, k f −

p k ² 2,w

α

,

p(x)

.

•

,

Gα = h(

)

.

•

,

,

G

h

,

.

直交多項式系

(1)

•

.

•

(r 1 (x), r 2 (x), . . .)

k 6 = l

(r _k , r _l ) _w = 0

,

w

(

).

上記の直交関数系の定義は

∀k,(r

, r

)

= 0

,

.

• Fourier

, { 1, sin x, cos x, sin 2x, cos 2x, . . . }

.

•

(r 1 (x), r 2 (x), . . .)

,

∀ k, (r k , r k ) w > 0

,

w

.

直交多項式系

(3)

• { 1, x, x ² , . . . }

Gram-Schmidt

(Legendre

, http://mathworld.wolfram.com/LegendrePolynomial.html ),

.

•

,

,

.

⊲ Legendre

• Legendre

(

)

,

p n (x)

.

p n (x) = 1 2 ⁿ n!

d ⁿ

dx ⁿ (x ² − 1) ⁿ , n = 0, 1, 2, . . . .

Legendre

, p 0 (x) = 1, p 1 (x) =

x, p 2 (x) = (3/2)x ² − (1/2)

.

直交多項式系

(5)

• Legendre

[ − 1, 1]

.

1

.

,

.

(p _k , p j ) = Z 1

− 1

p _k (x)p j (x)dx

•

, k 6 = j

(p k , p j ) = 0

.

•

,

,

, (p _k , p _k ) = 2/(2k + 1)

.

•

[

], pp. 157–158

.

,

I,

, 1980

.

直交多項式系

(7)

⊲ Laguerre

• Laguerre

(

)

,

l n (x)

.

l n (x) = e ^x n!

d ⁿ

dx ⁿ e ^−x x ⁿ

, n = 0, 1, 2, . . . .

• Leguerre

,

e ^−x

,

e ^x

,

.

• Leguerre

([

学辞典

]).

直交多項式系

(9)

• Leguerre

[0, ∞ )

.

e ^−x

.

,

.

(l _k , l j ) = Z _∞

0

l _k (x)l j (x)e ^−x dx

• k 6 = j

(l _k , l j ) = 0

, (l _k , l _k ) = (k!) ²

(

).

[

], pp. 278–279

.

• Fourier

Legendre

, Leguerre

[0, ∞ )

.

直交多項式系

(11)

⊲ Chebychev

• Chebychev

(

)

,

t n (x)

.

t _n (x) = cos(n arccos x), n = 0, 1, 2, . . . .

Chebychev

.

• n = 0

, t 0 (x) = cos 0 = 1

.

• n = 1

, t 1 (x) = cos(arccos x) = x

.

•

, k ≤ n

t k

x

k

, t _n +1 (x)

x

n + 1

.

直交多項式系

(13)

• arccos x = θ

,

cos(α + β) + cos(α − β ) = 2 cos α cos β

α = nθ, β = θ

, cos((n + 1)θ)) + cos((n − 1)θ) = 2 cos nθ cos θ

,

x = cos θ

, t n

, t n+1 (x)+ t _n− 1 (x) = 2xt n

となる

.

, t n+1 (x)

x

n + 1

.

• Chebychev

[ − 1, 1]

.

1/ √

1 − x ²

.

,

.

(t k , t j ) = Z 1

− 1

t k (x)t j (x) dx

√ 1 − x ² .

直交多項式系

(15)

• θ = arccos x (x = cos θ)

, t _k = cos kθ, t l = cos lθ, √

1 − x ² = sin θ, dx = − sin θdθ

,

(t _k , t _j ) = Z 0

π

cos kθ cos lθ − sin θdθ sin θ

となる

.

, Fourier

よって

, k 6 = j

(t _k , t _j ) = 0

, (t 0 , t 0 ) = π,

k ≥ 1

(t k , t k ) = π/2

.

• Chebychev

t k (x)

, x j = cos jπ

k , 0 ≤ j ≤

k

,

x j

t _k (x j )

1

,

t _k (x _j )

t _k (x _j +1 )

.

, Chebychev

.

直交多項式系

(17)

•

,

.

•

,

Fourier

.

• Scilab

Cheby- chev

.

chepol

Chebychev

.

•

chepol

, 2

(

,

)

,

Chebychev

.

•

.

-->chepol(3,’x’) ans =

3 - 3x + 4x -->chepol(5,’y’)

ans =

3 5

5y - 20y + 16y

• -->

Scilab

.

•

,

, Scilab

(-->)

.

Chebychev

• Chebychev

Lagrange

Chebychev

.

• Lagrange

あり

,

.

,

Chebychev

が知られている

([

,

]).

以下しばらくJ. Stoer, R. Bulirsch, Introduction to Numerical Analysis, 3/e,

Springer, 2002

R. Butt, Introduction to Numerical Analysis using MATLAB, Infinity Science Press, 2008

•

n

{ (x k , y k ) : 1 ≤ k ≤ n }

(

)

,

.

•

,

,

.

最小二乗法による曲線のあてはめ

(2)

•

,

α

g(x, α)

, P n

k =1 (y k − g(x k , α)) ²

α

.

(

)

.

.

•

g(x, α)

.

(

121 ∼ 126

)

.

•

,

.

,

g(x, α)

.

• α

,

α

.

最小二乗法による曲線のあてはめ

(4)

• f

,

min α

n

X

k =1

(y k − g(x k , α)) ²

である

.

,

.

• g(x)

,

g(x) = P m

k =1 α k r k (x)

, r _k (x)

α _k

,

.

• y = (y 1 , . . . , y _n ) ^T , r(x) = (r 1 (x), . . . , r _m (x)),

α = (α 1 , . . . , α m )

.

,

m 6 = n

.

最小二乗法による曲線のあてはめ

(6)

• R =





 r(x 1 )

.. . r(x _n )







. k y − Rα k ² 2

α

,

.

• k y − Rα k ² 2 = y ^T y − 2y ^T Rα+α ^T R ^T Rα

.

• ^∂k y − Rα k ² ₂

∂ α

α

.

• 0 ^T = − y ^T R+α ^T R ^T R,

y ^T R

, R ^T Rα = Ry

.

•

R ^T Rα = Ry

.

.

•

,

.

Scilab

(1)

• Scilab

reglin

(

125

).

•

,

, x

y

, y ≃ Ax + b

.

,

≃

.

• (x _i , y _i ) _i =1 ,...,N s

N s

(

), dim x = n x , dim y = n y

.

•

, A ∈ R ^{n y ×n x} , b ∈ R ^{n + y}

.

• reglin

, x i

N s

(

xmat

)

, y _i

N s

(

ymat

)

,

[A,b,s]=reglin[xmat,ymat]

.

• A, b

, s

.

Scilab

(3)

x

y

:

, x = 1, 2 . . . , 100, y = x + (

)

,

.

, a ≃ 1, b ≃ 0

.

x=1:100;

y=x+rand(1,100,’normal’);

[a,b,s]=reglin(x,y);

-->[a,b,s]

ans =

0.9997142 - 0.0073729 1.0491074

xmat

2

100

ymat

3

100

(正規分布)

xmat=rand(2,100,’normal’);

ymat=rand(3,100,’normal’);

[A,b,s]=reglin(xmat,ymat);

-->A A =

-0.0940391 -0.0510797 -0.0298481 0.0738067 -0.0712583 0.0198992

-->b b = -0.0986652 -0.0063689 -0.0370478

-->s s =

1.0819806 0.0205290 0.0777786

0.0205290 1.0073556 0.1069899

0.0777786 0.1069899 0.9121303

Scilab

(1)

• Scilab

小二乗法問題を解く関数

leastsq

,

.

•

,

y = f p (a, x) = a 0 + a 1 x + a 2 x ² + a 3 x ³

. x

− 1

1

0.1

, (a 0 , a 1 , a 2 , a 3 ) =

(1, 1, 1, 1)

, y

f p (a, x)

[ − 0.1, 0.1]

.

• leastsq

,

,

.

.

.

• leastsq

メータを探索するので

,

,

.

•

Scilab

.

Scilab

(3)

function y=fp(a,x) //パラメータを含む関数を定義 y=a(1)+a(2).*x+a(3).* (x.^2) + a(4).* (x.^3) endfunction

function e=fopt(a,xmat,ymat) //最小二乗法で使う関数 e=(norm(fp(a,xmat)-ymat))^2;

endfunction

aT=[1 2 3 4]; //パラメータの真値 xmat=-1:.1:1; //標本点の x

ymat=fp(aT,xmat)+0.1*(rand(xmat)-0.5); //y

a0=[1 1 1 1]; //パラメータ推定値の初期値

[f,aEst]=leastsq(list(fopt,xmat,ymat),a0); //最小化

実行結果の例:

->aEst aEst =

1.0050031 1.9672115 2.9836197 4.0484874

•

f

aEst

• leastsq

Scilab

•

x.^3

Scilab

(5)

•

datafit

.

leastsq

.

• leastsq

,

,

x

y

.

ページに

Scilab

(

じ問題

).

function y=fp(a,x) //パラメータを含む関数を定義,

y=a(1)+a(2).*x+a(3).*x.^2+a(4).*x.^3

endfunction

function e=fopt(a,z) //誤差関数, x

y

e=(z(2)-fp(a,z(1)))^2;

endfunction aT=[1;2;3;4];

xmat=-1:.1:1;

ymat=fp(aT,xmat)+0.1*rand(xmat,’normal’);

a0=[1;1;1;1]; //パラメータは列ベクトル

[aEst,err]=datafit(fopt,[xmat;ymat],a0); //返却値はパラメータ推定値と誤差

Scilab

(7)

•

, x

2

, y

, y ≃ a 0 + a ^T ₁ x + x ^T A 2 x

2

.

x

[0, 1]

.

•

Scilab

.

(fopt,xm,ym), a0);

(a0

function y=fp(a,x)

y=a(1)+a(2:3)’*x+x’*[a(4:5) a(6:7)]*x;

endfunction

function e=fopt(a,xm,ym) e=0;

for i=1:size(xm,2)

e=e+(ym(i)-fp(a,xm(:,i)))^2;

end endfunction aT=[1;1;1;1;0;0;1];

samples=100;

xm=rand(2,samples);

ym=[];

for i=1:samples ym=[ym fp(aT,xm(:,i))];

end

datafit

[aEst,err] = datafit (fopt, [xm;ym], a0);

(a0

function y=fp(a,x)

y=a(1)+a(2:3)’*x+x’*[a(4:5) a(6:7)]*x;

endfunction function e=fopt(a,z)

e=norm(z(3)-fp(a,z(1:2)))^2;

endfunction aT=[1;1;1;1;0;0;1];

samples=100;

xm=rand(2,samples);

ym=[];

for i=1:samples ym=[ym fp(aT,xm(:,i))];

end a0=rand(7,1);

[aEst,err]=datafit(fopt,[xm;ym],a0);

•

, x

y

2

, y ≃ a 0 + A 1 x + a 2 (x ^T x)

2

a 0

a ₂

2

, A 1

2

.

.

x

[0, 1]

.

•

Scilab

.

leastsq

[f,aEst] = leastsq (list (fopt,xm,ym), a0);

function y=fp(a,x)

y=a(1:2)+[a(3:4) a(5:6)]*x+a(7:8)*x’*x;

endfunction

function e=fopt(a,xm,ym) e=0;

for i=1:size(xm,2)

e=e+norm(ym(:,i)-fp(a,xm(:,i)))^2;

end endfunction aT=[1;1;1;0;0;1;1;1];

samples=100;

xm=rand(2,samples);

ym=[];

for i=1:samples ym=[ym fp(aT,xm(:,i))];

end

[xm;ym], a0);

function y=fp(a,x)

y=a(1:2)+[a(3:4) a(5:6)]*x+a(7:8)*x’*x;

endfunction function e=fopt(a,z)

e=norm(z(3:4)-fp(a,z(1:2)))^2;

endfunction aT=[1;1;1;0;0;1;1;1];

samples=100;

xm=rand(2,samples);

ym=[];

for i=1:samples ym=[ym fp(aT,xm(:,i))];

end a0=rand(8,1);

[aEst,err]=datafit(fopt,[xm;ym],a0);