4829-Gram-Schmidtの直交化法．.sagews

Gram-Schmidt orthogonalization

This is an example of Gram-Schmidt orthogonalization in a vector spaces defined by polynomials and integrals.

Sageを使って，線形代数（とくにGram-Schmidtの直交化法）のデモをした際のノートブックです．Sageをこれから始めたいWindowsユーザは、くろきげんさんの「Windows 7 ユーザーが数式処理ソフト Sage を使い始める方法」を参考にしてください．

 

まず変数$x$を定義する．

x=var('x')

内積を定義する．（念頭に置いているのは，たかだかn 次以下の実係数多項式の全体がなす実ベクトル空間）．

$$

\langle a, b \rangle := \int_0^1 a(x)b(x) \,dx.

$$

def innerprod_integration_0_1(a, b):
    return integrate(a*b, (x, 0, 1))

ノルムを定義する：

$

\|a \| := \sqrt{\langle a, a\rangle}.

$

def mynorm(a, inner_product_function):
    return sqrt(inner_product_function(a, a))

試しに、$1$, $x$, $x^2$のノルムを計算してみる．

mynorm(1, innerprod_integration_0_1)

$\renewcommand{\Bold}[1]{\mathbf{#1}}1$

mynorm(x, innerprod_integration_0_1)

$\renewcommand{\Bold}[1]{\mathbf{#1}}\sqrt{\frac{1}{3}}$

mynorm(x^2, innerprod_integration_0_1)

$\renewcommand{\Bold}[1]{\mathbf{#1}}\sqrt{\frac{1}{5}}$

Gram-Schmidtの直交化法．

 

内積空間$(V, \langle,\rangle)$を考える．一次独立なベクトルの組$(v_1,\dots, v_m)$が与えられたとき，次のようにしてorthonormal system $(e_1,\dots, e_m)$を定める手続きを，Gram-Schmidtの直交化法というのだった．

$

e_j := \frac{f_j}{\| f_j\|},\quad f_j = v_j - \sum_{k=1}^{j-1} \langle v_j, e_k\rangle e_k, \qquad (j = 1,\dots, m).

$

def gram_schmidt(vs, inner_product_function):
    ret = []
    for v in vs:
        f = v - sum([inner_product_function(v, r)*r for r in ret])
        ret.append(f/(mynorm(f, inner_product_function)))
    return(ret)

$P_3(\mathbf{R})$, たかだか3次の実係数多項式全体に，上で定義した内積で，Gram-Schmidtをやってみる．

rs = gram_schmidt([1,x,x^2, x^3], innerprod_integration_0_1)

for r in rs:
    print r.factor()

1 sqrt(3)*(2*x - 1) sqrt(5)*(6*x^2 - 6*x + 1) sqrt(7)*(10*x^2 - 10*x + 1)*(2*x - 1)

別の内積．積分区間を$[-1, 1]$にした．

def innerprod_integration_minus1_1(a, b):
    return integrate(a*b, (x, -1, 1))

いわゆるLegendreの多項式と呼ばれる，直交多項式の例．

legendre=gram_schmidt([1,x,x^2, x^3], innerprod_integration_minus1_1)

for r in legendre:
    print r.expand().simplify()

1/2*sqrt(2) 1/2*sqrt(3)*sqrt(2)*x 3/4*sqrt(5)*sqrt(2)*x^2 - 1/4*sqrt(5)*sqrt(2) 5/4*sqrt(7)*sqrt(2)*x^3 - 3/4*sqrt(7)*sqrt(2)*x

r

$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{7}{8} \, \sqrt{\frac{2}{7}} {\left(10 \, x^{3} - 3 \, \sqrt{3} \sqrt{2} \sqrt{\frac{2}{3}} x\right)}$

r.expand()

$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{35}{4} \, \sqrt{\frac{2}{7}} x^{3} - \frac{21}{8} \, \sqrt{3} \sqrt{2} \sqrt{\frac{2}{3}} \sqrt{\frac{2}{7}} x$

r.simplify()

$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{4} \, \sqrt{7} \sqrt{2} {\left(5 \, x^{3} - 3 \, x\right)}$

legendre_P(3,x)

$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{5}{2} \, {\left(x - 1\right)}^{3} + \frac{15}{2} \, {\left(x - 1\right)}^{2} + 6 \, x - 5$

_.simplify_full()

$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{5}{2} \, x^{3} - \frac{3}{2} \, x$

さらに別の内積の例．積分区間を$[-\pi, \pi]$にした．

def innerprod_integration_minus_pi_pi(a, b):
    return integrate(a*b, (x, -pi, pi))

教科書（S. Axler, Linear algebra done right, http://www.amazon.co.jp/dp/0387982582）の例（p.114-116）を参照．

es = gram_schmidt([1,x,x^2, x^3, x^4, x^5], innerprod_integration_minus_pi_pi)

for e in es:
    print e.factor()

1/2*sqrt(2)/sqrt(pi) 1/2*sqrt(3)*sqrt(2)*x/pi^(3/2) -1/4*sqrt(5)*sqrt(2)*(pi^2 - 3*x^2)/pi^(5/2) -1/4*sqrt(7)*sqrt(2)*(3*pi^2 - 5*x^2)*x/pi^(7/2) 3/16*sqrt(2)*(3*pi^4 - 30*pi^2*x^2 + 35*x^4)/pi^(9/2) 1/16*sqrt(11)*sqrt(2)*(15*pi^4 - 70*pi^2*x^2 + 63*x^4)*x/pi^(11/2)

Orthogonal projectionの計算

ベクトル空間$V$の部分空間$U$に対して，射影$P_U: V = U\oplus U^\bot \to U$をorthogonal projection（直交射影）というのだった．$f\in V$の直交射影$P_U(f)$は，$(e_1,\dots, e_m)$が$U$のorthonormal basisのとき，次のようにして計算できる：

$

P_U(f) := \sum_{j=1}^m \langle f, e_j\rangle e_j.

$

def orthogonal_projection(f, bases, inner_product_function):
    """This function computes the orthogonal projection of f onto the subspace spanned by bases."""
    return(sum([inner_product_function(f, base)*base for base in bases]))

$\sin(x)$の直交射影による多項式近似．

sin_approx=orthogonal_projection(sin(x), gram_schmidt([1,x,x^2,x^3,x^4,x^5], innerprod_integration_minus_pi_pi), innerprod_integration_minus_pi_pi)

sin_approx.expand().simplify_full()

$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{21 \, {\left(33 \, {\left(\pi^{4} - 105 \, \pi^{2} + 945\right)} x^{5} - 30 \, {\left(\pi^{6} - 125 \, \pi^{4} + 1155 \, \pi^{2}\right)} x^{3} + 5 \, {\left(\pi^{8} - 153 \, \pi^{6} + 1485 \, \pi^{4}\right)} x\right)}}{8 \, \pi^{10}}$

係数を近似値で見てみると次のようになる：

sum([RR(c[0])*x^c[1] for c in sin_approx.coefficients()])

$\renewcommand{\Bold}[1]{\mathbf{#1}}0.00564311797634682 \, x^{5} - 0.155271410633429 \, x^{3} + 0.987862135574674 \, x$

青線が多項式近似, 赤線が$\sin(x)$.　重なっていてほとんど誤差なし．

(plot(sin_approx, (x, -pi, pi))+plot(sin(x), (x, -pi, pi), color='red')).show()

テイラー展開との比較．

taylor(sin(x), x, 0, 5)

$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{120} \, x^{5} - \frac{1}{6} \, x^{3} + x$

直交射影による多項式近似（青線），sin(x) （赤線）, Taylor展開による近似（緑線）を重ねてプロットする．

(plot(sin_approx, (x, -pi, pi))+plot(sin(x), (x, -pi, pi), color='red')+plot(taylor(sin(x), x, 0, 5), (x, -pi, pi), color='green')).show()