数値解析
2023年度前期 第10週 [6月22日]

静岡大学大学院
工学研究科機械工学専攻
ロボット・計測情報分野
創造科学技術大学院情報科学専攻
三浦 憲二郎2次有理Bernstein基底(σ = 1/ 2) vs 1次一般化三角関数基底

講義アウトライン [6月15日]•連立1次方程式の反復解法•前回の復習•共役勾配法連立一次方程式の反復解法 p.93
•連立1次方程式
𝑨𝑨𝑨𝑨 = 𝒃𝒃•密(dense)行列と疎（sparse）行列
密行列：行列の要素のすべて，あるいはほとんどゼロでない行列
疎行列：行列の要素のほとんどがゼロである行列•疎行列•反復解法を用いるほうが速い．
•ヤコビ法とガウスザイデル法•共役勾配法

反復法の原理 p.93
連立1次方程式
𝑨𝑨𝑨𝑨 = 𝒃𝒃 (5.1)

を同値な方程式

𝒙𝒙 = 𝝓𝝓 𝒙𝒙 = 𝑴𝑴𝑴𝑴 + 𝑵𝑵𝑵𝑵に変形し，反復法により解を求める．𝑴𝑴を反復行列という．

縮小写像の原理

縮小写像の原理 定理5.1 p.94
𝑹𝑹𝒏𝒏 上で定義された関数𝒈𝒈(𝒙𝒙) が(1) 𝒙𝒙 ∈ 𝑹𝑹𝒏𝒏 → 𝒈𝒈(𝒙𝒙) ∈ 𝑹𝑹𝒏𝒏(2) 𝒙𝒙, 𝒚𝒚 ∈ 𝑹𝑹𝒏𝒏 → | 𝒈𝒈 𝒙𝒙 − 𝒈𝒈 𝒚𝒚 | ≤ 𝑳𝑳||𝒙𝒙 − 𝒚𝒚|| ∈ 𝑹𝑹𝒏𝒏
(3) 𝟎𝟎 ≤ 𝑳𝑳 < 𝟏𝟏を満たすとき，方程式
𝒙𝒙 = 𝒈𝒈(𝒙𝒙)の解𝒙𝒙 = 𝜶𝜶は𝑹𝑹𝒏𝒏 においてただ1つ存在する．そして，𝒙𝒙(𝟎𝟎) を初期値とする反復に𝒙𝒙(𝒌𝒌+𝟏𝟏) = 𝒈𝒈(𝒙𝒙(𝒌𝒌) )によって，生成される列{𝒙𝒙(𝒌𝒌) }は𝒌𝒌 → ∞のとき𝜶𝜶に収束する．

縮小写像の原理 定理5.2 p.94
あるノルム||・||について 𝐌𝐌 < 𝟏𝟏ならば反復法
𝒙𝒙(𝒌𝒌+𝟏𝟏) = 𝝓𝝓 𝒙𝒙 𝒌𝒌= 𝑴𝑴𝒙𝒙 𝒌𝒌 + 𝑵𝑵𝑵𝑵によって作られる列{𝒙𝒙(𝒌𝒌) }は(5.1)の解𝒙𝒙に収束する(証明)

任意の𝒙𝒙, 𝒚𝒚 ∈ 𝑹𝑹𝒏𝒏 に対して，
𝝓𝝓 𝒙𝒙 − 𝝓𝝓 𝒚𝒚= 𝑴𝑴𝑴𝑴 − 𝑴𝑴𝑴𝑴 ≤ 𝑴𝑴 � | 𝒙𝒙 − 𝒚𝒚 < | 𝒙𝒙 − 𝒚𝒚が成り立つので，𝝓𝝓 𝒙𝒙 は定理5.1の仮定を満たす．よって， {𝒙𝒙(𝒌𝒌) }は解𝒙𝒙に収束する．

行列Aの分解 p.95
𝑨𝑨 = 𝑫𝑫 + 𝑳𝑳 + 𝑼𝑼𝒂𝒂𝟏𝟏𝟏𝟏
𝑫𝑫 = ⋮𝟎𝟎
𝑳𝑳 = ⋮𝒂𝒂𝒏𝒏𝒏𝒏𝟎𝟎
𝑼𝑼 = ⋮𝟎𝟎⋯⋱⋮⋯ 𝒂𝒂𝒏𝒏𝒏𝒏
⋯ 𝟎𝟎
⋱ ⋮
⋯ 𝟎𝟎

⋯ 𝒂𝒂𝟏𝟏𝟏𝟏⋱⋮⋯ 𝟎𝟎

ヤコビ法とガウス・ザイデル法• ヤコビ法
とおいた反復法

𝑴𝑴 = −𝑫𝑫−𝟏𝟏 𝑳𝑳 + 𝑼𝑼 , 𝑵𝑵 = 𝑫𝑫−𝟏𝟏𝒙𝒙(𝒌𝒌+𝟏𝟏) = 𝑫𝑫−𝟏𝟏 {𝒃𝒃 − (𝑳𝑳 + 𝑼𝑼)𝒙𝒙(𝒌𝒌) }• ガウス・ザイデル法
とおいた反復法

𝑴𝑴 = − 𝑫𝑫 + 𝑳𝑳 −𝟏𝟏 𝑼𝑼, 𝑵𝑵 = 𝑫𝑫 + 𝑳𝑳 −𝟏𝟏𝒙𝒙(𝒌𝒌+𝟏𝟏) = 𝑫𝑫−𝟏𝟏 {𝒃𝒃 − 𝑳𝑳𝒙𝒙 𝒌𝒌+𝟏𝟏 − 𝑼𝑼𝒙𝒙 𝒌𝒌 }

共役勾配法 p.109
行列𝑨𝑨を対称正定値行列と仮定
正定値行列 p.57
n次正方行列𝑨𝑨が， 𝒙𝒙 ≠ 𝟎𝟎である任意のベクトル𝒙𝒙 ∈ 𝑹𝑹𝒏𝒏𝒙𝒙, 𝑨𝑨𝑨𝑨 > 𝟎𝟎.

• 𝑨𝑨が対称でない場合

𝑨𝑨𝒕𝒕 𝑨𝑨𝑨𝑨 = 𝑨𝑨𝒕𝒕 𝒃𝒃とすれば共役勾配法を使える．

共役勾配法 p.109
定理5.7 𝑨𝑨を𝒏𝒏次対称正定値行列とする．このとき，連立一次方程式の解�
𝒙𝒙は関数𝑨𝑨𝑨𝑨 = 𝒃𝒃

𝟏𝟏𝒇𝒇 𝒙𝒙 = 𝒙𝒙, 𝑨𝑨𝑨𝑨 − (𝒙𝒙, 𝒃𝒃)
𝟐𝟐を最小にする．ここで， �,� はベクトルの内積を表す．

共役勾配法 p.109
定理5.7の証明∀𝒉𝒉 ∈ 𝑹𝑹𝒏𝒏 𝒉𝒉 ≠ 𝟎𝟎 に対して
� + 𝒉𝒉 =𝒇𝒇 𝒙𝒙=𝟏𝟏� + 𝒉𝒉, 𝑨𝑨(�� + 𝒉𝒉, 𝒃𝒃𝒙𝒙
𝒙𝒙 + 𝒉𝒉) − 𝒙𝒙
𝟐𝟐𝟏𝟏𝟐𝟐�, 𝑨𝑨�𝒙𝒙𝒙𝒙) + (�𝒙𝒙, 𝑨𝑨𝑨𝑨) + (𝒉𝒉, 𝑨𝑨�𝒙𝒙) + (𝒉𝒉, 𝑨𝑨𝑨𝑨

� += 𝒇𝒇 𝒙𝒙

𝟏𝟏𝟐𝟐

�, 𝑨𝑨𝑨𝑨 + 𝒉𝒉, 𝑨𝑨�𝒙𝒙 + (𝒉𝒉, 𝑨𝑨𝑨𝑨) − 𝒉𝒉, 𝒃𝒃𝒙𝒙

�, 𝑨𝑨𝑨𝑨 = 𝑨𝑨�
𝒙𝒙, 𝒉𝒉 =ここで，𝑨𝑨は対称なので， 𝒙𝒙
よって，

� + 𝒉𝒉 = 𝒇𝒇 𝒙𝒙
� +𝒇𝒇 𝒙𝒙�, 𝒃𝒃 − 𝒉𝒉, 𝒃𝒃− 𝒙𝒙

𝟏𝟏𝟐𝟐𝒉𝒉, 𝑨𝑨�𝒙𝒙� +𝒉𝒉, 𝑨𝑨𝑨𝑨 + 𝒉𝒉, 𝑨𝑨�
𝒙𝒙 − 𝒉𝒉, 𝒃𝒃 = 𝒇𝒇 𝒙𝒙

𝑨𝑨の正定値性より 𝒉𝒉, 𝑨𝑨𝑨𝑨 > 𝟎𝟎なので，

𝟏𝟏𝟐𝟐

� + 𝒉𝒉 − 𝒇𝒇 𝒙𝒙
� > 𝟎𝟎, つまり 𝒇𝒇 𝒙𝒙� + 𝒉𝒉 > 𝒇𝒇 𝒙𝒙
�
𝒇𝒇 𝒙𝒙

𝒉𝒉, 𝑨𝑨𝑨𝑨共役勾配法 p.109,110
共役勾配法の原理
ベクトル 𝒑𝒑𝒊𝒊 と 𝒑𝒑𝒋𝒋 が𝒑𝒑𝒊𝒊 , 𝑨𝑨𝒑𝒑𝒋𝒋 = 𝟎𝟎,𝒊𝒊 ≠ 𝒋𝒋を満たすとき，ベクトル 𝒑𝒑𝒊𝒊 と 𝒑𝒑𝒋𝒋 は行列𝑨𝑨に関して共役である，あるいは𝑨𝑨に関して直交である，という．

𝒇𝒇 𝒙𝒙 を最小にする�
𝒙𝒙を求める𝒇𝒇(𝒙𝒙𝒌𝒌+𝟏𝟏 )は

𝜶𝜶𝒌𝒌 =

𝒙𝒙𝒌𝒌+𝟏𝟏 = 𝒙𝒙𝒌𝒌 + 𝜶𝜶𝒌𝒌 𝒑𝒑𝒌𝒌 ,(𝒌𝒌 = 𝟎𝟎, 𝟏𝟏, 𝟐𝟐, … )(𝒑𝒑𝒌𝒌 ,𝒓𝒓𝒌𝒌 )のとき最小となる．
(𝒑𝒑𝒌𝒌 ,𝑨𝑨𝒑𝒑𝒌𝒌 )ただし，誤差ベクトル𝒓𝒓𝒌𝒌+𝟏𝟏 = 𝒃𝒃 − 𝑨𝑨𝒙𝒙𝒌𝒌+𝟏𝟏 = 𝒃𝒃 − 𝑨𝑨 𝒙𝒙𝒌𝒌 + 𝜶𝜶𝒌𝒌 𝒑𝒑𝒌𝒌 = 𝒓𝒓𝒌𝒌 − 𝜶𝜶𝒌𝒌 𝑨𝑨𝒑𝒑𝒌𝒌

共役勾配法 p.110
探索方向𝒑𝒑𝒌𝒌+𝟏𝟏 はで定める．

𝒑𝒑𝒌𝒌+𝟏𝟏 = 𝒓𝒓𝒌𝒌+𝟏𝟏 + 𝜷𝜷𝒌𝒌 𝒑𝒑𝒌𝒌
(𝒓𝒓𝒌𝒌 , 𝑨𝑨𝒑𝒑𝒌𝒌 )𝜷𝜷𝒌𝒌 = −(𝒑𝒑𝒌𝒌 , 𝑨𝑨𝒑𝒑𝒌𝒌 )共役勾配法 Conjugate Gradient method: CG法
定理5.8 p.114 共役勾配法によって計算される誤差ベクトル𝒓𝒓𝒌𝒌 および方向ベクトル𝒑𝒑𝒌𝒌 について， 𝒊𝒊 ≠ 𝒋𝒋であれば，

𝒓𝒓𝒊𝒊 , 𝒓𝒓𝒋𝒋 = 𝟎𝟎, 𝒑𝒑𝒊𝒊 , 𝑨𝑨𝒑𝒑𝒋𝒋 = 𝟎𝟎.

定理5.9 p.115 𝑨𝑨は正則なn次対称正定値行列で，𝒃𝒃 はn次元ベクトルとする．
このとき，共役勾配法の反復をn回実行するすると𝑨𝑨𝑨𝑨 = 𝒃𝒃の解が得られる！

共役勾配法：プログラム その1/* 1ノルムの計算 */double vector_norm1(double* a, int m, int n){
int i;
double norm = 0.0;
for (i = m; i <= n; i++){norm += fabs(a[i]);}return norm;}/* Aベクトル a[m...n]とb[m...n]の内積を計算する */double inner_product(int m, int n, double* a, double* b){
int i;
double s = 0.0;
for (i = m; i <= n; i++) s += a[i] * b[i];return s;}/* 行列a[1...N][1...]とベクトルb[1...N]との積 c <- Ab */void matrix_vector_product(double** a, double* b, double* c){double wk;
int i, j;for (i = 1; i <= N; i++){wk = 0.0;
for (j = 1; j <= N; j++){wk += a[i][j] * b[j];}c[i] = wk;}}

共役勾配法：プログラム その2/* 共役勾配法 */double* cg(double** a, double* b, double* x){
double eps, * r, * p, * tmp, alpha, beta, work;
int i, k = 0;r = dvector(1, N); /* r[l...N] */p = dvector(1, N); /* p[l...N] */tmp = dvector(1, N); /* tmp[l...N] */matrix_vector_product(a, x, tmp); /* tmp <- A b */for( i = 1; i <= N; i++){p[i] = b[i] - tmp[i]; r[i] = p[i];}do{ /* alpha の計算 */matrix_vector_product(a, p, tmp); /* tmp <- A p_k */work = inner_product(1, N, p, tmp); /* work <- (p,Ap_k)*/alpha = inner_product(1, N, p, r) / work;/* x_{k+l} と r_{k+l}の計算 */for (i = 1; i <= N; i++) x[i] = x[i] + alpha * p[i];for (i = 1; i <= N; i++) r[i] = r[i] - alpha * tmp[i];/* 収束判定 */eps = vector_norm1(r, 1, N);
k++; /* 反復回数の更新 */if (eps < EPS) goto OUTPUT;/* beta と p_{k+l}の計算 */beta = -inner_product(1, N, r, tmp) / work;for (i = 1; i <= N; i++) p[i] = r[i] + beta * p[i];}while (k < KMAX);OUTPUT:;/* 領域の解放 */free_dvector(r, 1); free_dvector(p, 1); free_dvector(tmp, 1);
if (k == KMAX){
printf("答えが見つかりませんでした¥n");exit(1);}else {
printf("反復回数は%dです¥n", k); /* 反復回数を画面に表示 */return x;}}まとめ•連立1次方程式の反復解法•前回の復習•共役勾配法