広告 プログラミング

【解説＆実装】グラム・シュミットの正規直交化法

※本ページには、プロモーション（広告）が含まれています。

悩んでいる人
悩んでいる人

線形代数や数値解析の講義でグラム・シュミットの正規直交化法を学習したが、数値計算の実装に苦労している。

考え方や実装方法を教えて欲しい。

こんなお悩みを解決します。

今回は、グラム・シュミットの正規直交化法の解説とともに、C言語による実装を紹介していきたいと思います。

数式を交えつつ、考え方を順番に解説していきますので、興味がある方はぜひ最後までご覧ください。

グラム・シュミットの正規直交化法の手順

まずは、『筑波大学の武内修　先生のサイト』を参考に手順を整理します。また、簡単のため実ベクトルのみを対象とします。

前提

ここでは、与えられた一次独立な\(n\)個のベクトル\(\boldsymbol{a}_{1}, \boldsymbol{a}_{2}, \cdots, \boldsymbol{a}_{n}\)から正規直交系\(\boldsymbol{e}_{1}, \boldsymbol{e}_{2}, \cdots, \boldsymbol{e}_{n}\)を構築します。この時、\(\forall k (1 \le k \le n)\)について、\(\boldsymbol{e}_{k}\)は、\(\boldsymbol{a}_{1}, \boldsymbol{a}_{2}, \cdots, \boldsymbol{a}_{n}\)の一次結合で表せるものとします。

正規直交化の手順

まず、\(\boldsymbol{e}_{1}\)を以下のように設定します。

$$ \boldsymbol{e}_{1} = \frac{1}{\| \boldsymbol{a}_{1} \|}\boldsymbol{a}_{1} \quad (1) $$

次に、\(2 \le k \le n\)において、\(\boldsymbol{e}_{k}\)を以下のように求めます。

$$
\begin{eqnarray}
\boldsymbol{e}_{k} & = & \frac{1}{\| \boldsymbol{f}_{k} \|}\boldsymbol{f}_{k} \quad (2) \\
\text{ where } \boldsymbol{f}_{k} & = & \boldsymbol{a}_{k} - \sum_{i = 1}^{k - 1}(\boldsymbol{a}_{k}, \boldsymbol{e}_{i})\boldsymbol{e}_{i} \quad (3)
\end{eqnarray}
$$

ここで、\(\| \cdot \|\)はユークリッドノルム、\((\boldsymbol{u}, \boldsymbol{v}) = \boldsymbol{u}^{\top}\boldsymbol{v}\)です。

例題

一次独立な3個のベクトル\(\boldsymbol{a}_{1}, \boldsymbol{a}_{2}, \boldsymbol{a}_{3}\)が以下のように与えられているとし、これらのベクトルを正規直交化します。

$$
\begin{eqnarray*}
\boldsymbol{a}_{1} & = & (1, 1, 0)^{\top} \\
\boldsymbol{a}_{2} & = & (1, 0, 2)^{\top} \\
\boldsymbol{a}_{3} & = & (2, 1, 3)^{\top}
\end{eqnarray*}
$$

\(\boldsymbol{e}_{1}\)の計算

\(\| \boldsymbol{a}_{1} \| = \sqrt{2}\)ですので、式(1)より\(\boldsymbol{e}_{1}\)は、

$$
\begin{eqnarray*}
\boldsymbol{e}_{1} & = & \frac{1}{ \| \boldsymbol{a}_{1} \| }\boldsymbol{a}_{1} \\
& = & \frac{1}{ \sqrt{2} }(1, 1, 0)^{\top} \\
& = & \left( \frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}, 0 \right)^{\top}
\end{eqnarray*}
$$

となります。

\(\boldsymbol{e}_{2}\)の計算

\(k = 2\)の時の総和記号を展開すると、

$$
\boldsymbol{f}_{2} = \boldsymbol{a}_{2} - (\boldsymbol{a}_{2}, \boldsymbol{e}_{1})\boldsymbol{e}_{1}
$$

となります。ここで、

$$
\begin{eqnarray*}
(\boldsymbol{a}_{2}, \boldsymbol{e}_{1}) & = & (1, 0, 2)\left( \frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}, 0 \right)^{\top} \\
& = & \frac{1}{\sqrt{2}}
\end{eqnarray*}
$$

ですので、\(\boldsymbol{f}_{2}\)は、

$$
\begin{eqnarray*}
\boldsymbol{f}_{2} & = & \boldsymbol{a}_{2} - (\boldsymbol{a}_{2}, \boldsymbol{e}_{1})\boldsymbol{e}_{1} \\
& = & (1, 0, 2)^{\top} - \frac{1}{\sqrt{2}}\left( \frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}, 0 \right)^{\top} \\
& = & (1, 0, 2)^{\top} - \left( \frac{1}{2}, \frac{1}{2}, 0 \right)^{\top} \\
& = & \left( \frac{1}{2}, -\frac{1}{2}, 2 \right)^{\top}
\end{eqnarray*}
$$

となります。また、\(\| \boldsymbol{f}_{2} \| = \frac{3}{\sqrt{2}} \)です。よって、\(\boldsymbol{e}_{2}\)は、

$$
\begin{eqnarray*}
\boldsymbol{e}_{2} & = & \frac{1}{\| \boldsymbol{f}_{2} \|}\boldsymbol{f}_{2} \\
& = & \frac{\sqrt{2}}{3}\left( \frac{1}{2}, -\frac{1}{2}, 2 \right)^{\top} \\
& = & \left( \frac{\sqrt{2}}{6}, -\frac{\sqrt{2}}{6}, \frac{2\sqrt{2}}{3} \right)^{\top}
\end{eqnarray*}
$$

となります。

\(\boldsymbol{e}_{3}\)の計算

\(k = 3\)の時の総和記号を展開すると、

$$
\boldsymbol{f}_{3} = \boldsymbol{a}_{3} - (\boldsymbol{a}_{3}, \boldsymbol{e}_{1})\boldsymbol{e}_{1} - (\boldsymbol{a}_{3}, \boldsymbol{e}_{2})\boldsymbol{e}_{2}
$$

となります。ここで、

$$
\begin{eqnarray*}
(\boldsymbol{a}_{3}, \boldsymbol{e}_{1}) & = & (2, 1, 3)\left( \frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}, 0 \right)^{\top} \\
& = & \frac{3}{\sqrt{2}} \\
(\boldsymbol{a}_{3}, \boldsymbol{e}_{2}) & = & (2, 1, 3)\left( \frac{\sqrt{2}}{6}, -\frac{\sqrt{2}}{6}, \frac{2\sqrt{2}}{3} \right)^{\top} \\
& = & \frac{13\sqrt{2}}{6}
\end{eqnarray*}
$$

ですので、\(\boldsymbol{f}_{3}\)は、

$$
\begin{eqnarray*}
\boldsymbol{f}_{3} & = & \boldsymbol{a}_{3} - (\boldsymbol{a}_{3}, \boldsymbol{e}_{1})\boldsymbol{e}_{1} - (\boldsymbol{a}_{3}, \boldsymbol{e}_{2})\boldsymbol{e}_{2} \\
& = & (2, 1, 3)^{\top} - \frac{3}{\sqrt{2}}\left( \frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}, 0 \right)^{\top} - \frac{13\sqrt{2}}{6}\left( \frac{\sqrt{2}}{6}, -\frac{\sqrt{2}}{6}, \frac{2\sqrt{2}}{3} \right)^{\top} \\
& = & \left( 2 - \frac{3}{2} - \frac{13}{18}, 1 - \frac{3}{2} + \frac{13}{18}, 3 - \frac{26}{9} \right)^{\top} \\
& = & \left( -\frac{2}{9}, \frac{2}{9}, \frac{1}{9} \right)^{\top}
\end{eqnarray*}
$$

となります。また、\(\| \boldsymbol{f}_{3} \| = \frac{1}{3} \)です。よって、\(\boldsymbol{e}_{3}\)は、

$$
\begin{eqnarray*}
\boldsymbol{e}_{3} & = & \frac{1}{\| \boldsymbol{f}_{3} \|}\boldsymbol{f}_{3} \\
& = & 3\left( -\frac{2}{9}, \frac{2}{9}, \frac{1}{9} \right)^{\top} \\
& = & \left( -\frac{2}{3}, \frac{2}{3}, \frac{1}{3} \right)^{\top}
\end{eqnarray*}
$$

となります。

C言語による実装

実装結果を以下に格納しました。必要に応じて参照してください。

https://github.com/yuruto-free/gram-schmidt-orthonormalization

前提

引数として与えるベクトルのデータの構造は、以下のようになります。ここで、\(d\)は次元数、\(n\)はデータの数、\(\boldsymbol{a}_{k}=(a_{k, 1}, a_{k, 2}, \cdots, a_{k, d})^{\top}\)です。

$$
\begin{eqnarray*}
V = \left[
\begin{array}{cccc}
a_{1,1} & a_{1,2} & \cdots & a_{1,d}\\
a_{2,1} & a_{2,2} & \cdots & a_{2,d}\\
\vdots & \vdots & \ddots & \vdots \\
a_{n,1} & a_{n,2} & \cdots & a_{n,d}\\
\end{array}
\right]
\end{eqnarray*}
$$

すなわち、一次独立な\(n\)個のベクトル\(\boldsymbol{a}_{1}, \boldsymbol{a}_{2}, \cdots, \boldsymbol{a}_{n}\)を行ベクトル形式で並べたものが引数として与えるベクトルのデータ構造となります。

また、上記の変数名と関数の引数名の対応を以下に示します。

対象変数名関数の引数
次元数\(d\)dim
ベクトルの数\(n\)num
ベクトル\(V\)vecs
変数名と関数の引数の対応関係

実装のコア部分

グラム・シュミットの正規直交化法を実現しているのは、以下の部分となります。

int32_t GS_orthonormalization(int32_t dim, int32_t num, double *vecs) {
    int32_t ret = (int32_t)GS_NG;
    int32_t func_val;
    int32_t k;
    double norm;
    double *target;
    double *work = NULL;

    if (NULL != vecs) {
        work = (double *)malloc(sizeof(double) * (num - 1));
        if (NULL == work) {
            goto EXIT_GS_ORTH;
        }

        target = &vecs[0];
        /* e_{0}のノルムを計算 */
        func_val = euclid_norm(dim, target, &norm);
        if (((int32_t)RETURN_OK_GS != func_val) || (fabs(norm) < (double)MEPS_GS)) {
            goto EXIT_GS_ORTH;
        }
        /* e_{0}を計算 */
        func_val = scaling(dim, 1.0 / norm, target);
        if ((int32_t)RETURN_OK_GS != func_val) {
            goto EXIT_GS_ORTH;
        }
        /* e_{1}からe_{n-1}を計算 */
        for (k = 1; k < num; k++) {
            target = &vecs[k * dim];
            /* 直交ベクトルf_{k}を計算 */
            func_val = orthogonal_vector(dim, k, work, vecs);
            if ((int32_t)RETURN_OK_GS != func_val) {
                goto EXIT_GS_ORTH;
            }

            /* f_{k}のノルムを計算 */
            func_val = euclid_norm(dim, target, &norm);
            if (((int32_t)RETURN_OK_GS != func_val) ||
                (fabs(norm) < (double)MEPS_GS)) {
                goto EXIT_GS_ORTH;
            }
            /* e_{k}を計算 */
            func_val = scaling(dim, 1.0 / norm, target);
            if ((int32_t)RETURN_OK_GS != func_val) {
                goto EXIT_GS_ORTH;
            }
        }
        ret = (int32_t)GS_OK;
    }
EXIT_GS_ORTH:
    if (NULL != work) {
        free(work);
    }

    return ret;
}

// 中略

static int32_t orthogonal_vector(int32_t dim, int32_t idx, double *dots,
                                 double *vecs) {
    int32_t ret = (int32_t)RETURN_NG_GS;
    int32_t i;
    int32_t func_val;
    double *target;

    if ((NULL != dots) && (NULL != vecs)) {
        target = &vecs[idx * dim];

        /* 内積を計算 */
        for (i = 0; i < idx; i++) {
            func_val = inner_product(dim, target, &vecs[i * dim], &dots[i]);

            if ((int32_t)RETURN_OK_GS != func_val) {
                goto EXIT_ORTH_VEC;
            }
        }
        /* ベクトルの差を計算 */
        for (i = 0; i < idx; i++) {
            func_val = sub_vector(dim, dots[i], target, &vecs[i * dim]);

            if ((int32_t)RETURN_OK_GS != func_val) {
                goto EXIT_ORTH_VEC;
            }
        }
        ret = (int32_t)RETURN_OK_GS;
    }
EXIT_ORTH_VEC:

    return ret;
}

\(\boldsymbol{e}_{1}\)の計算

ベクトル\(V\)の1行目&vecs[0]を取り出し、ノルムを計算・正規化すると\(\boldsymbol{e}_{1}\)が求まります。ここで、C言語の配列は0から始まるので、対応関係に注意してください。

        target = &vecs[0];
        /* e_{0}のノルムを計算 */
        func_val = euclid_norm(dim, target, &norm);
        if (((int32_t)RETURN_OK_GS != func_val) || (fabs(norm) < (double)MEPS_GS)) {
            goto EXIT_GS_ORTH;
        }
        /* e_{0}を計算 */
        func_val = scaling(dim, 1.0 / norm, target);
        if ((int32_t)RETURN_OK_GS != func_val) {
            goto EXIT_GS_ORTH;
        }

\(\boldsymbol{e}_{k}\)の計算

\(\boldsymbol{e}_{2}\)以降は、一般的な形で計算します。実装方法は、以下のようになります。

  1. \(\boldsymbol{a}_{k}\)と\(\boldsymbol{e}_{1}, \boldsymbol{e}_{2}, \cdots, \boldsymbol{e}_{k - 1}\)のそれぞれの内積を計算します。
  2. 得られた内積を係数として、式(3)を計算します。
  3. 式(2)を用いて、得られた\(\boldsymbol{f}_{k}\)を正規化します。

1.のステップの実行時、配列vecsは以下のようになっています。

$$
\begin{eqnarray*}
\left[
\begin{array}{cccc}
e_{1,1} & e_{1,2} & \cdots & e_{1,d}\\
e_{2,1} & e_{2,2} & \cdots & e_{2,d}\\
\vdots & \vdots & \ddots & \vdots \\
e_{k - 1,1} & e_{k - 1,2} & \cdots & e_{k - 1,d}\\
a_{k,1} & a_{k,2} & \cdots & a_{k,d}\\
\end{array}
\right]
\end{eqnarray*}
$$

ここから、内積を計算する際は、作業用の配列を用意し、以下のように格納します。

$$
\begin{eqnarray*}
\boldsymbol{w} = \left[
\begin{array}{c}
(\boldsymbol{a}_{k}, \boldsymbol{e}_{1})\\
(\boldsymbol{a}_{k}, \boldsymbol{e}_{2})\\
\vdots \\
(\boldsymbol{a}_{k}, \boldsymbol{e}_{k - 1})\\
\end{array}
\right]
\end{eqnarray*}
$$

この結果を用いて、式(3)を計算します。ここで、上記の変数名と関数の引数の対応関係は以下のようになります。

対象変数名関数の引数名
正規直交化対象\(k\)idx
内積\(\boldsymbol{w}\)dots
ベクトル\(V\)vecs
変数名と関数の引数の対応関係
static int32_t orthogonal_vector(int32_t dim, int32_t idx, double *dots,
                                 double *vecs) {
    int32_t ret = (int32_t)RETURN_NG_GS;
    int32_t i;
    int32_t func_val;
    double *target;

    if ((NULL != dots) && (NULL != vecs)) {
        target = &vecs[idx * dim]; /* a_{k}に該当する変数 */

        /* 内積を計算 */
        for (i = 0; i < idx; i++) {
            /* 内積 (a_{k}, e_{i}) を計算 */
            func_val = inner_product(dim, target, &vecs[i * dim], &dots[i]);

            if ((int32_t)RETURN_OK_GS != func_val) {
                goto EXIT_ORTH_VEC;
            }
        }
        /* ベクトルの差を計算 */
        for (i = 0; i < idx; i++) {
            /* 式(3)の計算 */
            func_val = sub_vector(dim, dots[i], target, &vecs[i * dim]);

            if ((int32_t)RETURN_OK_GS != func_val) {
                goto EXIT_ORTH_VEC;
            }
        }
        ret = (int32_t)RETURN_OK_GS;
    }
EXIT_ORTH_VEC:

    return ret;
}

実行結果

main.cでは、以下のベクトルを正規直交化の対象としています。

$$
\begin{eqnarray*}
\left[
\begin{array}{ccc}
1 & 1 & 0\\
1 & 0 & 2\\
2 & 1 & 3\\
\end{array}
\right], \left[
\begin{array}{cccc}
1 & 1 & -2 & 2\\
0 & 1 & -1 & 0\\
3 & 5 & -2 & 1\\
\end{array}
\right], \left[
\begin{array}{ccccc}
1 & 2 & 1 & 3 & -1\\
0 & -2 & -3 & 3 & -2\\
2 & -3 & 0 & -4 & 3\\
3 & 1 & 2 & 3 & -4\\
\end{array}
\right]
\end{eqnarray*}
$$

テストコードを実行すると以下が得られます。少し誤差が含まれますが、おおむね計算できていることが分かります。

Dim: 3, Num: 3
[000] +0.70711 +0.70711 +0.00000
[001] +0.23570 -0.23570 +0.94281
[002] -0.66667 +0.66667 +0.33333
error: 4.1633363423443e-15

Dim: 4, Num: 3
[000] +0.31623 +0.31623 -0.63246 +0.63246
[001] -0.28604 +0.66742 -0.38139 -0.57208
[002] +0.67470 +0.51900 +0.51900 -0.07785
error: 3.6082248300318e-16

Dim: 5, Num: 4
[000] +0.25000 +0.50000 +0.25000 +0.75000 -0.25000
[001] -0.05000 -0.50000 -0.65000 +0.45000 -0.35000
[002] +0.85302 -0.36969 +0.06709 +0.05888 +0.35737
[003] +0.25084 -0.20149 +0.34284 -0.33565 -0.81626
error: 2.7384416402719e-03

スポンサードリンク



-プログラミング
-,