int gsl_linalg_QR_decomp (gsl_matrix *a, gsl_vector *tau)
を見てみます。
以下がソースですが、前回説明したように行列 gsl_matrix *a の各列ベクトルから求められる Householder 変換を *a 自身に施し続ける事で QR 分解を求めています。
01: #include <gsl/gsl_math.h>
02: #include <gsl/gsl_matrix.h>
03: #include <gsl/gsl_linalg.h>
04: #include <gsl/gsl_blas.h>
05:
06: int
07: gsl_linalg_QR_decomp (gsl_matrix *a, gsl_vector *tau)
08: {
09: const size_t m = a->size1;
10: const size_t n = a->size2;
11: if (tau->size != GSL_MIN (m, n)) GSL_ERROR ("size of tau must be MIN(m,n)", GSL_EBADLEN);
12: else {
13: size_t i;
14:
15: for (i = 0; i < GSL_MIN (m, n); i++) {
16: // a の i 列ベクトルの像
17: gsl_vector_view c_full = gsl_matrix_column (a, i);
18:
19: // a の i 列ベクトルの i 列以降の像
20: gsl_vector_view c = gsl_vector_subvector (&c_full.vector, i, m - i);
21:
22: /* c = a(i:m, i) を c から得られる Householder ベクトルで上書き。
23: 特に a(i, i) は beta = - sign (a(i, i)) * |c| で上書き。 */
24: double tau_i = gsl_linalg_householder_transform (&c.vector);
25:
26: gsl_vector_set (tau, i, tau_i);
27:
28: /* a の i 行 i + 1 列以降を Householder 変換。
29: a の i 列 については a(i, i) は beta、a(i + 1:m, i) = 0 となるが
30: a(i, i) = beta は 24 行の gsl_linalg_householder_transform で
31: 既に書き換え済み、a(i + 1:m, i) には Householder ベクトルの第2要素目が
32: 格納される(Householder ベクトルの第1要素は 1 で規定値なので省略)。 */
33: if (i + 1 < n) {
34: gsl_matrix_view b = gsl_matrix_submatrix (a, i, i + 1, m - i, n - (i + 1));
35: gsl_linalg_householder_hm (tau_i, &c.vector, &b.matrix);
36: }
37: }
38: }
39: return GSL_SUCCESS;
40: }
15 行から Houeholder 変換を入力された行列 gsl_matrix *a に適用して行きます。
i 回目の変換の前に 24 行で Householderベクトルを求めています。
この変換を i 列に施した場合、本来なら
となります。
a (i, i) がβになる事についてはここの最後の方で書いています。
本来上のようになるはずですが、Householderベクトルを求める際 i 列のベクトルの像を用いているので、元の行列の i 列 i 行以降の部分が計算された Householderベクトルで上書きされます。
なので i 列は飛ばして i + 1 列以降に i 回目の Householder 変換を施します。
となります。
この上三角部分が R となっています。
Q 自体は保存されませんが、各 Householder ベクトルの第 2 要素目以降（第 1 要素の値は 1 で既定値）が下三角部分に格納されていますので、これらから Q を再構成する事も可能で、以下の関数で Q 及び R が得られます。
int gsl_linalg_QR_unpack (const gsl_matrix *QR, const gsl_vector *tau, gsl_matrix *Q, gsl_matrix *R)
この関数では gsl_linalg_QR_decomp で分解済みの行列 gsl_matrix *QR と、同じ関数で返される τの値が格納された gsl_vector *tau を引数に取ります。
R は gsl_matrix *QR の上三角部分を取り出すだけです。
Q は gsl_matrix *QR の下三角部分に格納されたHouseholderベクトルと gsl_vector *tau の値を参照し gsl_linalg_householder_hm で単位行列に変換を施す事で得ています。
以下は Q と R を取り出すサンプルです。
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
void
fprintf_matrix (FILE *stream, gsl_matrix *a, const char *format)
{
int i, j;
for (i = 0; i < a->size1; i++) {
for (j = 0; j < a->size2; j++) {
fprintf (stream, format, gsl_matrix_get (a, i, j));
if (j < a->size2 - 1) fprintf (stream, "\t");
}
fprintf (stream, "\n");
}
return;
}
double data[] = {
1., 0., 2.,
0., 2., 0.,
2., 1., 1.
};
int
main (void)
{
gsl_matrix_view a = gsl_matrix_view_array (data, 3, 3);
gsl_vector *tau = gsl_vector_alloc (GSL_MIN (a.matrix.size1, a.matrix.size2));
gsl_matrix *q = gsl_matrix_alloc (a.matrix.size1, a.matrix.size1);
gsl_matrix *r = gsl_matrix_alloc (a.matrix.size1, a.matrix.size2);
/* QR分解 */
gsl_linalg_QR_decomp (&a.matrix, tau);
/* 行列 Q、R を構成 */
gsl_linalg_QR_unpack (&a.matrix, tau, q, r);
fprintf (stdout, "QR = [\n");
fprintf_matrix (stdout, &a.matrix, "%.3f");
fprintf (stdout, "]\n\nQ = [\n");
fprintf_matrix (stdout, q, "%.3f");
fprintf (stdout, "]\n\nR = [\n");
fprintf_matrix (stdout, r, "%.3f");
fprintf (stdout, "]\n");
gsl_vector_free (tau);
gsl_matrix_free (q);
gsl_matrix_free (r);
return 0;
}
この出力の結果は
QR = [
-2.236 -0.894 -1.789
0.000 -2.049 0.293
0.618 0.110 -1.309
]
Q = [
-0.447 0.195 -0.873
0.000 -0.976 -0.218
-0.894 -0.098 0.436
]
R = [
-2.236 -0.894 -1.789
0.000 -2.049 0.293
0.000 0.000 -1.309
]
となりました。
まぁ、直接 Q や R を取り出す必要はあまりないとは思いますが。