修正Cholesky分解

[] [] []


修正Cholesky分解のアルゴリズム

A=RDR

変数を節約するためrii=diとする。さらにrji=sij(i<j)と使ってもよい。

【外積形式】

/* Aの上三角をRに入れる */

for i=1 to n

  for j=1 to i-1

        rij = 0

    for j=i to n

        rij = aij

 

/* Rの計算 */

for k=1 to n

    for i=k+1 to n

        for j=i to n

            rij = rij - rki* rkj/rkk

    for j=k+1 to n

        rkj = rkj/rkk

【内積形式行ごと】

/* Aの上三角をRに入れる */

for i=1 to n

    for j=1 to i-1

        rij = 0

    rii = aii

    for j=i+1 to n

        sij = aij

/* Rの計算 */

for k=1 to n

    for i=1 to k-1

        rkk = rkk - rik* sik

    for j=k+1 to n

        for i=1 to k-1

            skj = skj - rik* sij

        rkj = skj/rkk

【内積形式列ごと】

/* Aの上三角をRに入れる */

for i=1 to n

    for j=1 to i-1

        rij = 0

    rii = aii

    for j=i+1 to n

        sij = aij

/* Rの計算 */

for k=1 to n

    for i=1 to k-1

        for j=1 to i-1

            sik = sik - rji* sjk

        rik = sik/rii

    for i=1 to k-1

        rkk = rkk - rik* sik

連立1次方程式Ax=bを解くには

/* Ry=byについて解く(前進消去) */

for i=1 to n

   tmp = 0

   for j=1 to i-1

      tmp = tmp + rji* yj

   yi = bi - tmp

/* Rx=D-1yxについて解く(後退代入) */

xn = yn / rnn

for i=n-1 to 1

   tmp = 0

   for j=i+1 to n

      tmp = tmp + rji* xj

   xi = yi / rii - tmp

正定値のエルミート行列Aを三角行列に分解するアルゴリズムを導く。

A対角要素aiiはエルミート行列だから実数であり

さらに正定値だからaii>0であることに注意。

まずnn列のA11列目とそれ以外に分けて次のように書く

対角要素a11の下を消去すると

さらにa11の右側も消去できて、A3つの行列に分解できる。

ここで

このn-1n-1列の行列も正定値のエルミート行列になっている。なぜなら

の右辺は明らかに正定値のエルミート行列だからである。

したがって同様に分解を続ければ、対角行列D対角要素が1である下三角Rと上三角Rに分解される。これを修正コレスキー分解という

A=RDR

具体的には、次の手順になる。ただし変数を節約するためrii=diとする。

手順

コード

Aの上三角を作業領域Rに入れる。

・ステップk (k=1 to n)において

 k+1行目から下のrを求める

  rij = rij - rki* rkj/rkk (i=k+1 to n, j=i to n)

  (LU分解がj=k+1からnまで計算するのに比べ半分で済む)

 k行目のrを求める

  rkj = rkj/rkk (i=k+1 to n)

【外積形式】

/* Aの上三角をRに入れる */

for i=1 to n

    for j=1 to i-1

        rij = 0

    for j=i to n

        rij = aij

 

/* Rの計算 */

for k=1 to n

    for i=k+1 to n

        for j=i to n

            rij = rij - rki* rkj/rkk

    for j=k+1 to n

        rkj = rkj/rkk

◇◇

上記のアルゴリズムは外積aaを使うので外積形式とよばれるが、内積形式のアルゴリズムもある

アルゴリズムを導く方針は、A=RDRと分解できたとして要素ごとに式を書き下し、これをd,rについて解く。

簡単のためn=3の場合を書き下すと

d,rについて解く手順として、行ごとか列ごとかの2とおりがある。

◇◇◇

行ごとに解くと

1行目 d1=a11 , r12=a12/d1 , r13=a13/d1

2行目 d2=a22 - d1 r12*r12 , r23=(a23 - d1 r12*r13)/d2

3行目 d3=a33 - d1 r 13*r13 - d2 r23*r23

この手順を一般化すると次のようになる。ただし変数を節約するためrii=diとする。

手順

コード

Aの上三角を作業領域Rに入れる。

・ステップk (k=1 to n)において

k行目のrkkを求める

k行目のrkjを求める

ここでriiで割ってから再びかけるのは無駄なので

中間変数sij=riirijを使うとよい

【内積形式,行ごと】

/* Aの上三角をRに入れる */

for i=1 to n

    for j=1 to i-1

        rij = 0

    rii = aii

    for j=i+1 to n

        sij = aij

/* Rの計算 */

for k=1 to n

    for i=1 to k-1

        rkk = rkk - rik* sik

    for j=k+1 to n

        for i=1 to k-1

            skj = skj - rik* sij

        rkj = skj/rkk

◇◇◇◇

列ごとに解くと

1列目 d1=a11

2列目 r12=a12/d1 ,d2=a22 - d1 r12*r12

3列目 r13=a13/d1 , r23=(a23 - d1 r12*r13)/d2 , d3=a33 - d1 r 13*r13 - d2 r23*r23

この手順を一般化すると次のようになる。ただし変数を節約するためrii=diとする。

手順

コード

Aの上三角を作業領域Rに入れる。

・ステップk (k=1 to n)において

k列目のrikを求める

k列目のrkkを求める

ここでriiで割ってから再びかけるのは無駄なので

中間変数sij=riirijを使うとよい

【内積形式,列ごと】

/* Aの上三角をRに入れる */

for i=1 to n

    for j=1 to i-1

        rij = 0

    rii = aii

    for j=i+1 to n

        sij = aij

/* Rの計算 */

for k=1 to n

    for i=1 to k-1

        for j=1 to i-1

            sik = sik - rji* sjk

        rik = sik/rii

    for i=1 to k-1

        rkk = rkk - rik* sik


[] [] []