今回は2重対角化と、正方対称行列の3重対角化について。
これらの変換は
int gsl_linalg_bidiag_decomp (gsl_matrix *a, gsl_vector *tau_u, gsl_vector *tau_v)
及び
int gsl_linalg_symmtd_decomp (gsl_matrix *a, gsl_vector *tau)
で求められます。
2重対角化とは、行列を対角成分と上対角(super diagonal)成分のみが非ゼロとなるようにするものです。
この変換は特異値分解の前処理などに使われます。
3重対角化とは、相似変換で対角成分と上・下(sub)対角成分のみ非ゼロにするものです。
これは固有値計算の際に用いられたりします。
これらの変換も Householder変換で求められます。
まず2重対角化ですが、まず行列の1列目についてHouseholder変換を施します。
こうすると1列目の [ \beta_11, 0, ..., 0]^t の部分を変更することなしに1行3列目以降を0に出来ます。
この操作を繰り返すことで行列を2重対角化できます。
以下が gsl_linalg_bidiag_decomp のソースです。
01: #include <gsl/gsl_math.h>
02: #include <gsl/gsl_matrix.h>
03: #include <gsl/gsl_blas.h>
04:
05: /*
06: 行列 gsl_matrix *a の2重対角化を計算
07: a->size1 >= a->size2 を仮定
08: 列についてのHouseholder変換係数 (size2 個) が gsl_vector *tau_u に
09: 行についてのHouseholder変換係数 (size2 - 1 個) が gsl_vector *tau_v に格納される。
10: */
11: int
12: cl_linalg_bidiag_decomp (gsl_matrix *a, gsl_vector *tau_u, gsl_vector *tau_v)
13: {
14: if (a->size1 < a->size2) GSL_ERROR ("bidiagonal decomposition requires M>=N", GSL_EBADLEN);
15: if (tau_u->size != a->size2) GSL_ERROR ("size of tau_U must be N", GSL_EBADLEN);
16: if (tau_v->size + 1 != a->size2) GSL_ERROR ("size of tau_V must be (N - 1)", GSL_EBADLEN);
17: {
18: size_t i;
19: const size_t m = a->size1;
20: const size_t n = a->size2;
21:
22: for (i = 0; i < n; i++) {
23:
24: // i 列ベクトルによるHouseholder変換
25: {
26: gsl_vector_view c = gsl_matrix_column (a, i);
27: // i 列ベクトルの i 行以降の像
28: gsl_vector_view v = gsl_vector_subvector (&c.vector, i, m - i);
29: // i 列ベクトルの i 行以降をHouseholderベクトルで上書き
30: double tau_i = gsl_linalg_householder_transform (&v.vector);
31:
32: // i 行 i + 1 列以降の部分行列をHouseholder変換
33: if (i + 1 < n) {
34: gsl_matrix_view s = gsl_matrix_submatrix (a, i, i + 1, m - i, n - (i + 1));
35: gsl_linalg_householder_hm (tau_i, &v.vector, &s.matrix);
36: }
37:
38: // tau_u の第 i 成分に変換係数を格納
39: gsl_vector_set (tau_u, i, tau_i);
40:
41: }
42:
43: // i 行ベクトルによるHouseholder変換
44: if (i + 1 < n) {
45: gsl_vector_view r = gsl_matrix_row (a, i);
46: // i 行ベクトルの i + 1 列以降の像
47: gsl_vector_view v = gsl_vector_subvector (&r.vector, i + 1, n - (i + 1));
48: // i 行ベクトルの i + 1 列以降をHouseholderベクトルで上書き
49: double tau_i = gsl_linalg_householder_transform (&v.vector);
50:
51: // i + 1 行 i + 1 列以降の部分行列をHouseholder変換
52: if (i + 1 < m) {
53: gsl_matrix_view s = gsl_matrix_submatrix (a, i + 1, i + 1, m - (i + 1), n - (i + 1));
54: gsl_linalg_householder_hm (tau_i, &v.vector, &s.matrix);
55: }
56:
57: // tau_v の第 i 成分に変換係数を格納
58: gsl_vector_set (tau_v, i, tau_i);
59:
60: }
61: }
62:
63: }
64: return GSL_SUCCESS;
65: }
14行目で行列サイズのチェックが行われていますが、この関数では 行数 (a->size1) <= 列数 (->size2) であるものしか扱えません。
そうでない場合は転置して下2重対角行列を求めることになるようです。
この関数内でやっていることは、まさに上に書いたことそのままです。
この関数で上書きされた gsl_matrix *a の対角成分、上対角成分が 2重対角化後の対角、上対角成分に、その他がこの変換のHouseholderベクトルになります。
この関数で返された gsl_matrix *a、gsl_vector *tau_u, gsl_vector *tau_v を用い
int gsl_linalg_bidiag_unpack2 (gsl_matrix *a, gsl_vector *tau_u, gsl_vector *tau_v, gsl_matrix *p)
で、2重対角化後の対角、上対角成分が tau_u, tau_v に格納されます。
また、列についてのHouseholder変換行列を Q、行についてのHouseholder変換行列を P とすると
Q' * A * P
で2重対角化されますが、Q が a を上書き、P が p に格納されます。
tau_u、tau_v は a から値を抜き出す事で、Q、P は a に格納されているHouseholderベクトルによるHouseholder変換を単位行列に施すことで求めています。
ちなみに・・・
有名な LAPACK では、特異値分解を計算する dgesvd_ などでは、2重対角化された行列を求める以下の関数
void dgebrd_ (int *m, int *n, double *a, int *lda, double *d, double *e, double *tauq, double *taup,
double *work, int *lwork, int *info)
を呼び出しています。
int *m, *n が行列 double *a の行・列数
int *lda がリーディングディメンジョン（こいつの効能については次回話すとしましょう・・・）
double *d, *e には対角成分、上(m < n なら下)対角成分が格納
double *tauq, *taup は gsl の場合の tau_u, tau_v に当たります。
gsl は row major ですが、LAPACK では col major で行列を格納することに注意して・・・
double *a に行列を格納して
int m = a の行数
int n = a の列数
int lda = m; (とりあえず）
int min_mn = MIN (m, n);
double *d = (double *) malloc (min_mn * sizeof (double));
double *e = (double *) malloc ((min_mn - 1) * sizeof (double));
double *tauq = (double *) malloc (min_mn * sizeof (double));
double *taup = (double *) malloc (min_mn * sizeof (double));
double wkopt;
double *work;
lwork = -1;
dgebrd_ (&m, &n, a, &lda, d, e, tauq, taup, &wkopt, &lwork, &info);
if (info != 0) エラー;
lwork = (int) wkopt;
work = (double *) malloc (lwork * sizeof (double));
dgebrd_ (&m, &n, a, &lda, d, e, tauq, taup, work, &lwork, &info);
これで gsl の関数と同じ結果が得られます。
長くなったので3重対角化は次回で。。。