ヤコビ(Jacobi)法

対称行列A(n×n)のn個の固有値・固有ベクトルを求める方法.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
/*!
 * Jacobi法による固有値の算出
 * @param[inout] a 実対称行列.計算後,対角要素に固有値が入る
 * @param[out] v 固有ベクトル(aと同じサイズ)
 * @param[in] n 行列のサイズ(n×n)
 * @param[in] eps 収束誤差
 * @param[in] iter_max 最大反復回数
 * @return 反復回数
 */
int EigenJacobiMethod(float *a, float *v, int n, float eps = 1e-8, int iter_max = 100)
{
    float *bim, *bjm;
    float bii, bij, bjj, bji;
 
    bim = new float[n];
    bjm = new float[n];
 
    for(int i = 0; i < n; ++i){
        for(int j = 0; j < n; ++j){
            v[i*n+j] = (i == j) ? 1.0 : 0.0;
        }
    }
 
    int cnt = 0;
    for(;;){
        int i, j;
 
        float x = 0.0;
        for(int ia = 0; ia < n; ++ia){
            for(int ja = 0; ja < n; ++ja){
                int idx = ia*n+ja;
                if(ia != ja && fabs(a[idx]) > x){
                    i = ia;
                    j = ja;
                    x = fabs(a[idx]);
                }
            }
        }
 
        float aii = a[i*n+i];
        float ajj = a[j*n+j];
        float aij = a[i*n+j];
 
        float alpha, beta;
        alpha = (aii-ajj)/2.0;
        beta  = sqrt(alpha*alpha+aij*aij);
 
        float st, ct;
        ct = sqrt((1.0+fabs(alpha)/beta)/2.0);    // sinθ
        st = (((aii-ajj) >= 0.0) ? 1.0 : -1.0)*aij/(2.0*beta*ct);    // cosθ
 
        // A = PAPの計算
        for(int m = 0; m < n; ++m){
            if(m == i || m == j) continue;
 
            float aim = a[i*n+m];
            float ajm = a[j*n+m];
 
            bim[m] =  aim*ct+ajm*st;
            bjm[m] = -aim*st+ajm*ct;
        }
 
        bii = aii*ct*ct+2.0*aij*ct*st+ajj*st*st;
        bij = 0.0;
 
        bjj = aii*st*st-2.0*aij*ct*st+ajj*ct*ct;
        bji = 0.0;
 
        for(int m = 0; m < n; ++m){
            a[i*n+m] = a[m*n+i] = bim[m];
            a[j*n+m] = a[m*n+j] = bjm[m];
        }
        a[i*n+i] = bii;
        a[i*n+j] = bij;
        a[j*n+j] = bjj;
        a[j*n+i] = bji;
 
        // V = PVの計算
        for(int m = 0; m < n; ++m){
            float vmi = v[m*n+i];
            float vmj = v[m*n+j];
 
            bim[m] =  vmi*ct+vmj*st;
            bjm[m] = -vmi*st+vmj*ct;
        }
        for(int m = 0; m < n; ++m){
            v[m*n+i] = bim[m];
            v[m*n+j] = bjm[m];
        }
 
        float e = 0.0;
        for(int ja = 0; ja < n; ++ja){
            for(int ia = 0; ia < n; ++ia){
                if(ia != ja){
                    e += fabs(a[ja*n+ia]);
                }
            }
        }
        if(e < eps) break;
 
        cnt++;
        if(cnt > iter_max) break;
    }
 
    delete [] bim;
    delete [] bjm;
 
    return cnt;
}

ハウスホルダー変換

実対称行列を三重対角行列に変換する.さらに三重対角行列から固有値/固有ベクトルを求めるまでをハウスホルダー法と呼ぶ. ハウスホルダー変換はQR法におけるヘッセンベルグ行列への変換にも用いられる.

ハウスホルダー変換では固有値を求めたい行列Aから対称な直交行列Pを定義し, 1行1列目が三重対角行列になった行列PAPを計算,これを繰り返す. Pが対称な直交行列(P=P^-1=P^T)なのでAとPAPの固有値は等しい.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
/*!
 * ハウスホルダー変換で実対称行列を三重対角行列に変換
 * @param[inout] a 元の行列(n×n).変換された三重対角行列を格納
 * @param[in] n 行列のサイズ
 */
void HouseholderTransformation(float *a, int n)
{
    float *u, *v, *q;
    u = new float[n];
    v = new float[n];
    q = new float[n];
 
    for(int j = 0; j < n; ++j){
        u[j] = 0.0;
        v[j] = 0.0;
        q[j] = 0.0;
    }
 
    for(int k = 0; k < n-2; ++k){
        // sの計算
        float s = 0;
        for(int i = k+1; i < n; ++i){
            s += a[i*n+k]*a[i*n+k];
        }
        s = -((a[(k+1)*n+k] >= 0) ? 1 : -1)*sqrt(s);
 
        // |x-y|の計算
        float alpha = sqrt(2.0*s*(s-a[(k+1)*n+k]));
        if(fabs(alpha) < 1e-8) continue;
 
        // uの計算
        u[k+1] = (a[(k+1)*n+k]-s)/alpha;
        for(int i = k+2; i < n; ++i){
            u[i] = a[i*n+k]/alpha;
        }
 
        // Auの計算
        q[k] = alpha/2.0;
        for(int i = k+1; i < n; ++i){
            q[i] = 0.0;
            for(int j = k+1; j < n; ++j){
                q[i] += a[i*n+j]*u[j];
            }
        }
 
        // v=2(Au-uu^T(Au))の計算
        alpha = 0.0;
        for(int i = k+1; i < n; ++i){
            alpha += u[i]*q[i];
        }
        v[k] = 2.0*q[k];
        for(int i = k+1; i < n; ++i){
            v[i] = 2.0*(q[i]-alpha*u[i]);
        }
 
        // A = PAP = A-uv^T-vu^Tの計算
        a[k*n+(k+1)] = a[(k+1)*n+k] = s;
        for(int i = k+2; i < n; ++i){
            a[k*n+i] = 0.0;
            a[i*n+k] = 0.0;
        }
        for(int i = k+1; i < n; ++i){
            a[i*n+i] = a[i*n+i]-2*u[i]*v[i];
            for(int j = i+1; j < n; ++j){
                a[i*n+j] = a[i*n+j]-u[i]*v[j]-v[i]*u[j];
                a[j*n+i] = a[i*n+j];
            }
        }
    }
 
    delete [] u;
    delete [] v;
    delete [] q;
}

QR法

正方行列A(n×n)の固有値と固有ベクトルを求める.Aは複素行列でも可.

QR法ではまず,非対称行列をハウスホルダー変換によりヘッセンベルグ行列に変換する (対称行列ではないので三重対角行列ではないことに注意). ハウスホルダー変換によるヘッセンベルグ行列への変換は以下.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
/*!
 * ハウスホルダー変換で行列をヘッセンベルグ行列に変換
 * @param[inout] a 元の行列(n×n).変換されたヘッセンベルグ行列を格納
 * @param[in] n 行列のサイズ
 */
void HouseholderTransformationForQR(float *a, int n)
{
    float *b, *p, *q;
    float *u;
    b = new float[n*n];
    p = new float[n*n];
    q = new float[n*n];
    u = new float[n];
 
    for(int k = 0; k < n-2; ++k){
        // sの計算
        float s = 0;
        for(int i = k+1; i < n; ++i){
            s += a[i*n+k]*a[i*n+k];
        }
        s = -((a[(k+1)*n+k] >= 0) ? 1 : -1)*sqrt(s);
 
        // |x-y|の計算
        float alpha = sqrt(2.0*s*(s-a[(k+1)*n+k]));
        if(fabs(alpha) < 1e-8) continue;
 
        // uの計算
        u[k+1] = (a[(k+1)*n+k]-s)/alpha;
        for(int i = k+2; i < n; ++i){
            u[i] = a[i*n+k]/alpha;
        }
 
        // Pの計算
        for(int i = k+1; i < n; ++i){
            for(int j = i; j < n; ++j){
                if(j == i){
                    p[i*n+j] = 1.0-2.0*u[i]*u[i];
                }
                else{
                    p[i*n+j] = -2.0*u[i]*u[j];
                    p[j*n+i] = p[i*n+j];
                }
            }
        }
 
        // PAの計算
        for(int i = k+1; i < n; ++i){
            for(int j = k+1; j < n; ++j){
                q[i*n+j] = 0.0;
                for(int m = k+1; m < n; ++m){
                    q[i*n+j] += p[i*n+m]*a[m*n+j];
                }
            }
        }
 
        // A = PAP^Tの計算
        for(int i = 0; i <= k; ++i){
            b[i*n+k] = a[i*n+k];
        }
        b[(k+1)*n+k] = s;
        for(int i = k+2; i < n; ++i){
            b[i*n+k] = 0.0f;
        }
        for(int j = k+1; j < n; ++j){
            for(int i = 0; i <= k; ++i){
                b[i*n+j] = 0.0;
                for(int m = k+1; m < n; ++m){
                    b[i*n+j] += a[i*n+m]*p[j*n+m];
                }
            }
            for(int i = k+1; i < n; ++i){
                b[i*n+j] = 0.0;
                for(int m = k+1; m < n; ++m){
                    b[i*n+j] += q[i*n+m]*p[j*n+m];
                }
            }
        }
 
        for(int j = 0; j < n*n; ++j){
            a[j] = b[j];
        }
    }
 
    delete [] u;
    delete [] b;
    delete [] p;
    delete [] q;
}

行列Aをヘッセンベルグ行列Hへ変換後,QR法による固有値を算出する. ハウスホルダー変換後の行列を引数hに渡す.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
/*!
 * QR法による固有値の算出
 * @param[inout] h 元の行列(ヘッセンベルグ行列).計算後,対角要素に固有値が入る
 * @param[in] n 行列のサイズ(n×n)
 * @param[in] eps 収束誤差
 * @param[in] iter_max 最大反復回数
 * @return 反復回数
 */
void QR(float *h, int n, float eps = 1e-8, int iter_max = 100)
{
    float *r, *q, *t;
    r = new float[n*n];
    q = new float[n*n];
    t = new float[n*n];
 
    float *u, *v;
    u = new float[n];
    v = new float[n];
 
    // R = H : ヘッセンベルグ行列
    for(int i = 0; i < n; ++i){
        for(int j = 0; j < n; ++j){
            int idx = i*n+j;
            r[idx] = h[idx];
        }
    }
 
    int cnt = 0;
    for(;;){
        // Q=I (単位行列)
        for(int i = 0; i < n; ++i){
            for(int j = 0; j < n; ++j){
                q[i*n+j] = (i == j) ? 1.0 : 0.0;
            }
        }
 
        for(int k = 0; k < n-1; ++k){
            // sinθ,cosθの計算
            float alpha = sqrt(r[k*n+k]*r[k*n+k]+r[(k+1)*n+k]*r[(k+1)*n+k]);
            if(fabs(alpha) < 1e-8) continue;
 
            float c = r[k*n+k]/alpha;
            float s = -r[(k+1)*n+k]/alpha;
 
            // Rの計算
            for(int j = k+1; j < n; ++j){
                u[j] = c*r[k*n+j]-s*r[(k+1)*n+j];
                v[j] = s*r[k*n+j]+c*r[(k+1)*n+j];
            }
            r[k*n+k] = alpha;
            r[(k+1)*n+k] = 0.0;
            for(int j = k+1; j < n; ++j){
                r[k*n+j]     = u[j];
                r[(k+1)*n+j] = v[j];
            }
 
            // Qの計算
            for(int j = 0; j <= k; ++j){
                u[j] = c*q[k*n+j];
                v[j] = s*q[k*n+j];
            }
            q[k*n+(k+1)] = -s;
            q[(k+1)*n+(k+1)] = c;
            for(int j = 0; j <= k; ++j){
                q[k*n+j]     = u[j];
                q[(k+1)*n+j] = v[j];
            }
        }
 
        // RQの計算
        for(int i = 0; i < n; ++i){
            for(int j = 0; j < n; ++j){
                float rq = 0.0;
                for(int m = 0; m < n; ++m){
                    rq += r[i*n+m]*q[j*n+m];
                }
 
                t[i*n+j] = rq;
            }
        }
 
        // 収束判定
        float e = 0.0;
        for(int i = 0; i < n; ++i){
            e += fabs(t[i*n+i]-h[i*n+i]);
        }
        if(e < eps) break;
        if(cnt > iter_max) break;
 
        for(int i = 0; i < n; ++i){
            for(int j = 0; j < n; ++j){
                int idx = i*n+j;
                r[idx] = t[idx];
                h[idx] = t[idx];
            }
        }
 
        cnt++;
    }
 
    delete [] r;
    delete [] q;
    delete [] t;
    delete [] u;
    delete [] v;
}

QR関数を使用した例を示す. 行列

 3  0  0
-2 -2  4
 0 -1  3

の固有値を求める.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
    int n = 3;
    float *h = new float[n*n];
    float h0[3][3] = { {3, 0, 0}, 
                       {-2, -2, 4}, 
                       {0, -1, 3} };
 
    for(int i = 0; i < n; ++i){
        for(int j = 0; j < n; ++j){
            h[i*n+j] = h0[i][j];
        }
    }
 
    HouseholderTransformationForQR(h, n);
    QR(h, n, 1e-8, 100);
 
    cout << "QR" << endl;
    for(int i = 0; i < n; ++i){
        for(int j = 0; j < n; ++j){
            cout << h[i*n+j] << " ";
        }
        cout << endl;
    }
 
    cout << "eigen values : " << h[0] << ", " << h[4] << ", " << h[8] << endl;
    
    delete [] h;

結果は,

     |3 -0.333333 -0.29814|
QR : |3.66271e-007 2 -5.36656|
     |0 6.50788e-011 -1|
eigen values : 3, 2, -1

数学的の求めてみる.特性方程式は,

(λ-3)((λ+2)(λ-3)+4) = 0

変形すると

(λ-3)(λ-2)(λ+1) = 0

となり,λ=3,2,-1で上記結果と一致する.


トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2012-06-06 (水) 18:25:23 (3396d)