[C言語]グラムシュミットの直交化のQR分解を用いて連立方程式を解く[コード付]

どうも　シュモクザメです。
今日はQR分解について解説してみたいと思います。
まあ行列の変換ですね。
以前にガウスの消去法やLU分解も解説しているので是非どうぞ！

QR分解とは？

まずQR分解っていうのは、LU分解と同様に連立方程式を解くために係数行列を変形しようってノリですね。

係数行列ってなんや

係数行列っていうのはその名前の通りですね。
以下のように連立方程式を行列の計算に変換することを行列形式表示といいます。

一般的にするとこんな感じ

Aが係数が入っている係数行列、
ｂが定数項が入ってる定数項ベクトル、
ｘは求める解が入っている解ベクトルです。

このAをA=QRに分解しようっていうのがQR分解です。
このときにQ,Rを
Qは直交行列
Rは上三角行列になるようにします。

ここまできたら元の計算はこのようになります。

このRX=Q^tbまで変形できれば後退代入法で解くことができます！

後退代入法、、、？

聞きなれないですよね。
後退代入法っていうのは、係数行列が上三角行列の時にすぐにXを計算できる方法です。
こんな式

この時に解が以下のようになります。

なんでこうなるのかは分かりませんが、せっかくあるんで使わせてもらいましょう。
この形ならプログラミングで一瞬で解けますね。

QR分解のメリット

うーん、QR分解で変形すれば後退代入法を使って解けるのはわかったけど、別にガウスの消去法でよくない？？

いい質問ですね〜
確かに手間を考えるとどっちも変わらんやんってなりますよね。

行列計算の際には条件数condという概念があります。
簡単に説明すると誤差の大きさのことです。
条件数が大きいほど誤差が大きいってことなのですが、ガウスの消去法とQR分解の条件数を比べるとQR分解の方が条件数が小さいんです。
つまり誤差が少ないってことなので一般的には良条件ってわけですね。
まあ条件数の計算方法とかに関してはまたググってみてください

一般的なQとR

QとRの各要素を求める前に、ここでは係数行列Aと分解行列Q,Rのi列j行の要素を記事中ではそれぞれa_ij、q_ij、r_ijとします。（図中ではijは添字です）

QとRは以下のように求められます。

これはグラムシュミットの直交化という方法を使って求めているのですが、
この計算に関しても特に解説しません。
（そういうものだと思ってありがたく使わせてもらいましょう）

コードの実装

まず色々な関数を作って準備しましょう

準備1

行列の表示の関数です。

void SHOW_1(double B[N]){//1次元配列の中身を表示する
    int i;
    for(i=0; i<N; i++){
        printf("%f\n",B[i]);
    }
}

void SHOW_2(double A[N][N]){//2次元配列の中身を表示する
    printf("///////\n");
    int i, j;
    for(i=0; i<N; i++){
        for(j=0; j<N; j++){
            printf("%f ",A[i][j]);
        }
        printf("\n");
    }
}

まあそのまんまですね。
次は転置と積の関数です。

void transpose(double A[N][N],double B[N][N]){//Aの転置をBに代入
    
    for(int i=0; i<N; i++){
        for(int j=0; j<N; j++)
            B[j][i] = A[i][j];
    }
}

void multi_NxN_Nx1(double A[N][N], double B[N], double C[N]){//AとBの積をCに代入
    for(int i=0; i<N; i++){
        for(int j=0; j<N; j++){
            C[i] += A[i][j]*B[j];
        }
    }
}

これも特に解説しなくても大丈夫そうですね。

準備2

次はベクトル単位で操作する際の関数です。
上の画像で示した通り、QとRを求める際はAやQを列ごとにベクトルとしたものを多用します。
それをコード上で表現するとa_kはA[k][i](0<=i<N)となることに注意しましょう。

内積を計算する関数

まずは内積の計算です。

double inner_product(double A[N][N], int a, double E[N][N], int e){//Aのa列目とEのe列目の内積を求める
    double sum=0;
    for(int i=0; i<N; i++){
        sum += A[i][a]*E[i][e];
    }
    
    return sum;
}

内積の定義が分かっていれば簡単だと思います。
内積自体は整数ですが、途中の計算で小数が出てくるのでdouble型で返しましょう

ベクトルの長さを求める関数

次はベクトルの長さです。

double norm(double A[N][N], int a){//Aのa列目の長さを求める
    double sum=0;
    for(int i=0; i<N; i++)
        sum += A[i][a]*A[i][a];
    
    return sqrt(sum);
}

sqrtを用いるのでmath.hをインクルードしときましょう。

あるベクトルにスカラー倍した他のベクトルを足す関数

void plus(double A[N][N], int a, double B[N][N], int b, double T){//AにT倍したBを代入
    for(int i=0; i<N; i++)
        A[i][a] = A[i][a] + T*B[i][b];
}

見出しの通りの機能を持つ関数です。めちゃくちゃ重要な関数です。

あるベクトルをスカラー倍する関数

void scalar_times(double A[N][N], int a, double T){//Aのa列目をT倍する
    for(int i=0; i<N; i++)
        A[i][a] = T*A[i][a];
}

これも見出しの通りの機能を持つ関数d(ry

グラムシュミットの直交化を用いてQを求める関数。

ながーい準備がやっと終わりました。ここから本番ですね。

void GS_ortho(double A[N][N], double Q[N][N], double Q_t[N][N]){//グラムシュミットの直交化
    
    for(int k=0; k<N; k++){
        plus(Q,k,A,k,1);
        for(int i=0; i<k; i++){
            plus(Q, k, Q_t, i, -inner_product(A, k, Q_t, i));
            printf("inner %f",-inner_product(A, k, Q_t, i));
        }
        plus(Q_t, k, Q, k, 1/norm(Q,k));
    }
}

少しややこしいのでじっくりやっていきましょう。

仮引数

仮引数は
Aが係数行列が入っている配列
Qはこれから直交化を計算したものを代入する配列
Q_tはこれからQを正規化したものを代入する配列です

最終的に使うのはQ_tだから仮引数もAとQ_tだけで良いんじゃいの？

と思うんですが、めんどくさいことに直交化しただけの行列がRの計算で必要なんですね。
てことで橋渡し役として仮引数にQも入れてあげてます。

外側のfor

外側のforはkを0からN-1まで動かします。
kは上の画像で説明した、q_kを計算するよってことです。
おさらいですがq_kはQのk列目のベクトルですね。

その後関数plusでQ_kにA_kを代入してます。
この部分ですね。

このときスカラー倍のところは1倍としてます。

内側のfor

内側のforはiを0からk-1まで動かして、第二項以降を計算してます。
この部分ですね。

計算式通りに関数plusの引数を与えてます。
スカラー倍のところを-innner_productとして内積を与えてるのがポイントです。

最後の関数plus

ここは直交化が完了したベクトル(Q[i][k])を自身の長さで割って正規化させている部分ですね(Q_t[i][k])
スカラー倍を1/normとしてます。

これでQが求まりました！！

Rを求める関数

次はRです。と言っても上で示した通りに要素を計算するだけです。
引数はA,Q,Q_tは上と全く同じで、RはそのままRを代入する用の配列です。

void cal_R(double A[N][N], double Q[N][N],double Q_t[N][N], double R[N][N]){//Rの計算
    for(int i=0; i<N; i++)
        R[i][i] = norm(Q, i);

    for(int i=N-2; i>=0; i--){
        for(int j=i+1; j<N; j++){
            R[i][j] = inner_product(A, j, Q_t, i);
        }
    }
}

一つ目のfor

一つ目のforはRの対角成分を計算してます。

Q、つまり正規化するまえのQの長さをnormを用いて代入してます。

2つ目のfor

ここは残りの部分を計算します。

流れとしては一番下の要素(N-1行目)から初めて、1行ずつ上にずらして求めて行ってます。
外側のforはiをN-2から0まで動かします。
内側のforはjをi+1からN-1まで動かします。
ちょっとややこしいですが以下のイメージで走査してますね。

内側の処理では関数innner_productを用いて内積を計算してR[i][j]に代入してます。

これでRも求まりましたね！

後退代入法

後退代入法に関しては以前ガウスの消去法について書いた時に用いたものを流用します。
と言っても特に難しいところはありません。

void ANS(double A[N][N], double B[N], double X[N]){//後退代入法
    int i,j;
    X[N-1] = B[N-1]/A[N-1][N-1];
    
    for(i=N-2; i>=0; i--){
        
        double S=0;
        for(j=i+1; j<N; j++)
            S = S + A[i][j] * X[j];
        
        X[i] = (B[i]-S)/A[i][i];
    }
}

Aが係数行列、Bが定数項ベクトル、Xが解ベクトルです。
Aには上で求めたRが、
Bには上で求めたQの転置と、元の問題のbの積のベクトルを入れます。

全体のコード

#include <stdio.h>
#include <math.h>

#define N 3

void transpose(double A[N][N],double B[N][N]){//Aの転置をBに代入
    
    for(int i=0; i<N; i++){
        for(int j=0; j<N; j++)
            B[j][i] = A[i][j];
    }
}

void multi_NxN_Nx1(double A[N][N], double B[N], double C[N]){//AとBの積をCに代入
    for(int i=0; i<N; i++){
        for(int j=0; j<N; j++){
            C[i] += A[i][j]*B[j];
        }
    }
}

void SHOW_1(double B[N]){//1次元配列の中身を表示する
    int i;
    for(i=0; i<N; i++){
        printf("%f\n",B[i]);
    }
}

void SHOW_2(double A[N][N]){//2次元配列の中身を表示する
    printf("///////\n");
    int i, j;
    for(i=0; i<N; i++){
        for(j=0; j<N; j++){
            printf("%f ",A[i][j]);
        }
        printf("\n");
    }
}

double inner_product(double A[N][N], int a, double E[N][N], int e){//Aのa列目とEのe列目の内積を求める
    double sum=0;
    for(int i=0; i<N; i++){
        sum += A[i][a]*E[i][e];
    }
    
    return sum;
}

void plus(double A[N][N], int a, double B[N][N], int b, double T){//AにT倍したBを代入
    for(int i=0; i<N; i++)
        A[i][a] = A[i][a] + T*B[i][b];
}

double norm(double A[N][N], int a){//Aのa列目の長さを求める
    double sum=0;
    for(int i=0; i<N; i++)
        sum += A[i][a]*A[i][a];
    
    return sqrt(sum);
}

void scalar_times(double A[N][N], int a, double T){//Aのa列目をT倍する
    for(int i=0; i<N; i++)
        A[i][a] = T*A[i][a];
}

void GS_ortho(double A[N][N], double Q[N][N], double Q_t[N][N]){//グラムシュミットの直交化
    
    for(int k=0; k<N; k++){
        plus(Q,k,A,k,1);
        for(int i=0; i<k; i++){
            plus(Q, k, Q_t, i, -inner_product(A, k, Q_t, i));
            printf("inner %f",-inner_product(A, k, Q_t, i));
        }
        plus(Q_t, k, Q, k, 1/norm(Q,k));
    }
}

void cal_R(double A[N][N], double Q[N][N],double Q_t[N][N], double R[N][N]){//Rの計算
    for(int i=0; i<N; i++)
        R[i][i] = norm(Q, i);

    for(int i=N-2; i>=0; i--){
        for(int j=i+1; j<N; j++){
            R[i][j] = inner_product(A, j, Q_t, i);
            printf("%d %d :%f\n",i,j,R[i][j]);
        }
    }
}

void ANS(double A[N][N], double B[N], double X[N]){//後退代入法
    int i,j;
    X[N-1] = B[N-1]/A[N-1][N-1];
    
    for(i=N-2; i>=0; i--){
        
        double S=0;
        for(j=i+1; j<N; j++)
            S = S + A[i][j] * X[j];
        
        X[i] = (B[i]-S)/A[i][i];
    }
}

int main(int argc, const char * argv[]) {
    double A[N][N] = { {2,3,4},
                    {1,5,3},
                    {3,0,2}
                };
    
    double B[N] = {2,1,4};
    
    SHOW_2(A);
    
    double Q[N][N]={0};
    double Q_n[N][N]={0};
    double R[N][N]={0};
    
    GS_ortho(A, Q, Q_n);
    printf("Qn\n");
    SHOW_2(Q_n);
    printf("Q\n");
    SHOW_2(Q);
    
    cal_R(A, Q, Q_n, R);
    printf("R\n");
    SHOW_2(R);
    
    double Q_t[N][N]={0};
    transpose(Q_n, Q_t);
    printf("Qt\n");
    SHOW_2(Q_t);
    
    double B_t[N]={0};
    multi_NxN_Nx1(Q_t, B, B_t);
    
    double X[N]={0};
    ANS(R, B_t, X);
    
    printf("X\n");
    SHOW_1(X);
    
    return 0;
}