を、行列 A の左から掛けたものを求める関数は
int gsl_linalg_householder_hm (double tau, gsl_vector *v, gsl_matrix *a)
となります。
ここで H の i 行 j 列成分は
なので A' = H・A の i 行 j 列成分は
式（１）
となります。
ここで i = 1 の時
なので、式（１）は式（２）
となります。
さらに i = 1 か否かで
式（３）
式（４）
となります。
以下が gsl_linalg_householder_hm のソースです。
gsl/linalg/householder.c
001: #include <gsl/gsl_math.h>
002: #include <gsl/gsl_matrix.h>
003: #include <gsl/gsl_blas.h>
004: #include <gsl/gsl_linalg.h>
005:
006: int
007: gsl_linalg_householder_hm (double tau, gsl_vector *v, gsl_matrix *a)
008: {
009: if (tau == 0.) return GSL_SUCCESS;
010:
011: {
012: // Householderベクトルの第2要素以下
013: gsl_vector_const_view v1 = gsl_vector_const_subvector (v, 1, v->size - 1);
014: // 入力行列の2行目以降
015: gsl_matrix_view a1 = gsl_matrix_submatrix (a, 1, 0, a->size1 - 1, a->size2);
016: size_t j;
017:
018: for (j = 0; j < a->size2; j++) {
019: double wj = 0.;
020: // 行列 a の j列ベクトルの第2要素以下
021: gsl_vector_view a1j = gsl_matrix_column (&a1.matrix, j);
022: double a0j = gsl_matrix_get (a, 0, j);
023:
024: /* 式（２）の wj */
025: // wj = sum_{k=1}^{n-1} v[k] * a[k, j] : 式（２）のsummation
026: gsl_blas_ddot (&v1.vector, &a1j.vector, &wj);
027: // wj = a[0, j] + sum_{k=1}^{n-1} v[k] * a[k, j] : 式（２）
028: wj += a0j;
029:
030: /* 式（３）: a[0, j] = a[0, j] - tau * wj */
031: gsl_matrix_set (a, 0, j, a0j - tau * wj);
032:
033: /* 式（４）: a[i, j] = a[i, j] - tau * wj * v[i] (i > 0) */
034: gsl_blas_daxpy (- tau * wj, &v1.vector, &a1j.vector);
035: }
036: return GSL_SUCCESS;
037: }
038: }
本来のコードでは BLAS を使用しない場合のものもあるのですが、BLASを使う部分だけ抜粋しています。
例えば、この関数を以下の様に呼び出した場合の動作を見てみます。
sample_main.c
001: #include <gsl/gsl_math.h>
002: #include <gsl/gsl_matrix.h>
003: #include <gsl/gsl_linalg.h>
004:
005: double data[] = {
006: 1., 0., 1.,
007: 0., 1., 1.,
008: 1., 1., 0. };
009:
010: /* 行列出力用の関数 */
011: void
012: fprintf_matrix (FILE *stream, gsl_matrix *a, const char *format)
013: {
014: int i, j;
015: const size_t m = a->size1;
016: const size_t n = a->size2;
017:
018: for (i = 0; i < m; i++) {
019: for (j = 0; j < n; j++) {
020: fprintf (stream, format, gsl_matrix_get (a, i, j));
021: if (j < n - 1) fprintf (stream, "\t");
022: }
023: fprintf (stream, "\n");
024: }
025: return;
026: }
027:
028: int
029: main (void)
030: {
031: gsl_matrix_view a = gsl_matrix_view_array (data, 3, 3);
032: gsl_vector_view view = gsl_matrix_column (&a.matrix, 0);
033: gsl_vector *v = gsl_vector_alloc (view.vector.size);
034: double tau;
035:
036: gsl_vector_memcpy (v, &view.vector);
037: tau = gsl_linalg_householder_transform (v);
038: gsl_linalg_householder_hm (tau, v, &a.matrix);
039:
040: fprintf_matrix (stdout, &a.matrix, "%.5f");
041:
042: return 0;
043: }このサンプルでは行列
を入力としています。
これは 数理ファイナンス様 で例として挙げられていた行列と同じものです。
ここでは上の行列 A の 1 列ベクトルから求められる Householder 変換を行列 A 自身に適用する事を行っています。
sample_main.c の 31 行目で上の行列の像を作成。
32 行でその 1 列目の列ベクトル （1, 0, 1)^t の像を作成し 36 行でそれを新たなベクトル v にコピー。
37 行の gsl_linalg_householder_transform で v からHouseholder ベクトルと τ を求める。
ここでは
τ = 1.707
v = [ -1.414, 0.0 0.414 ]^t
と、前回と同じ結果が得られます。
38行で上の τ、 v によるHouseholer変換を行列 A に左から施します。
gsl/linalg/householder.c の gsl_linalg_householder_hm に入って・・・
13 行目で Householder ベクトル v の第2要素以下のベクトルの像を作成します。
今の場合
v1.vector.data[] = { 0.0, 0.414 }
15 行目で入力行列の1行目以降の部分行列の像を作成
a1.matrix.data[] = {
0, 1., 1.,
1., 1., 0.
}
18 行目以降で入力行列の各列を変換後の要素で上書きしてゆきます。
j = 0 の場合。
21 行で a1 の j = 0 列ベクトルの像を作成。
a1j.vector.data[] = { 0., 1. }^t
26、27 行で式（２）の wj = a[0, 0] + v1.vector * a1j.vector を計算
wj = 1 + (0 + 0.414 * 1) = 1.414
31 行で式（３）に従い a[0, 0] → a[0, 0] - tau * wj= 1. - 1.7071 * 1.4142 = - 1.414 に更新
続く 34 行で式（４）に従い a[1,0]、 a[2,0] を更新
a[1, 0] → a[1, 0] - tau * wj * v[1] = 0 - tau * wj * 0 = 0
a[2, 0] → a[2, 0] - tau * wj * v[2] = 1 - 1.7071 * 1.4142 * 0.4142 = 0
で、35 行の j （= 0） のループを抜けたところで
となります。
これは v （= a の1列ベクトル） を鏡像変換で第一軸上のベクトルへ変換する操作に当たりますね。
ちなみに 34 行の axpy (a, x, y) は、a をスカラー、x、y をベクトルとして、ベクトル y = y + a * x を計算するものになります。
34 行の引数は a = - tau * wj、 x = v1.vector と y = a1j.vector ですので
a1j.vector → a1j.vector - tau * wj * v1.vector
と、a1j.vector の各要素を上書きしているのですが、 a1j.vector は入力行列 a の各要素（j = 0 の場合 1列の第2要素以下）をポイントしているので、 a1j.vector の上書きで a の各要素自体が書き換えられます。
この後 18 行に戻り j = 1, 2 について計算し、最終的に行列 a は
に置き換えられます。
これは 数理ファイナンス様 の途中計算の結果
に一致します。
なんか結構簡単な計算を馬鹿丁寧に追い過ぎた感もありますが・・・（汗
ところで、今回の様に行列の列ベクトルから得られる Householder 変換をその行列自身に適用する場合、例えば上の例の 1 列ベクトルのHouseholder変換は
となるので νは
となります。
この場合式（３）で j = 1 の場合、つまり 1 行 1 列目の要素は
となります。
j 列ベクトル（の j 行目以降）から得られる Householder変換を元の行に適用すると、変換された行列の j 行 j 列要素は βの値に一致します（ j 列 j + 1 行以降は 0 ）。
gsl_householder_transform では Householderベクトルの第一要素にβを代入（本来の値は1 ）していましたが、これはHouseholder変換で QR 分解を計算する際に効いてきます。