Cholesky分解

[←]　[↑]　[→]

Cholesky分解のアルゴリズム

A=R†R

 

【外積形式】

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

for i=1 to n

  for j=1 to n

    if(i<=j) rij = aij

    if(i>j )  rij = 0

/* Rの計算 */

for k=1 to n

   for i=k+1 to n

      for j=i to n

        rij = rij - rki* rkj/rkk

   rkk= sqrt(rkk)

   for j=k+1 to n

      rkj = rkj/rkk

【内積形式、行ごと】

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

for i=1 to n

  for j=1 to n

    if(i<=j) rij = aij

    if(i>j )  rij = 0

/* Rの計算 */

for k=1 to n

   for i=1 to k-1

      rkk = rkk - rik* rik

   rkk = sqrt(rkk)

   for j=k+1 to n

      for i=1 to k-1

         rkj = rkj - rik* rij

      rkj = rkj/rkk

【内積形式、列ごと】

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

for i=1 to n

  for j=1 to n

    if(i<=j) rij = aij

    if(i>j )  rij = 0

/* Rの計算 */

for k=1 to n

   for i=1 to k-1

      for j=1 to i-1

         rik = rik - rji* rjk

      rik = rik/rii

      rkk = rkk - rik* rik

   rkk = sqrt(rkk)

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

/* R†y=bをyについて解く(前進消去) */

for i=1 to n

   tmp = 0

   for j=1 to i-1

      tmp = tmp + rji* yj

   yi = (bi - tmp)/rii

/* Rx=yをxについて解く(後退代入) */

for i=n to 1

   tmp = 0

   for j=i+1 to n

      tmp = tmp + rij xj

   xi = (yi - tmp)/rii

前項の修正Cholesky分解の式(2)において、

対角要素a11の平方根をとって下三角と上三角にかけると

と対角行列は単位行列になるから、下三角と上三角の2つに分解される。

A=R†R

これがオリジナルのCholesky分解であるが、修正Cholesky分解と比べると平方根演算が増え、R†y=bの計算でも割り算が増えており、計算量で見ると損である。

[←]　[↑]　[→]