JavaScript
Node.js
ECMAScript

Nodeの為の不完全コレスキー分解付共役勾配ソルバ

More than 1 year has passed since last update.

はじめに

Nodeユーザーの為の共役勾配法のテンプレート。

収束しないこともあるらしい。
不完全コレスキー分解付共役勾配法は誤差の影響を減らしこの問題を解決する為にある。

圧力のポアソン方程式にディリクレ境界条件とノイマン境界条件を付加するように、
実装するときには計算過程のいたるところに条件分岐が入ることになる。

プログラムは下の記事と大学水工学研究室のものを読んでいる。

jicfus.jp

slis.tsukuba.ac.jp

共役勾配法

反復法として利用され、コレスキー分解のような直接法では大きすぎて取り扱えない、大規模な疎行列を解くために利用される。
そのような問題は偏微分方程式などを数値的に解く際に常に現れる。
= ja.wikipedia.org

xを解とするn次方程式を考える。

{\bf A} {\bf x} = {\bf b}

ここで、Aは正定値対称行列なので注意。

Node
const A = new Float64Array([
  1, 2, 4, 6,
  2, 3, 2, 8,
  4, 2, 2, 4,
  6, 8, 4, 1
])

const x = new Float64Array([
  0, 0, 0, 0
])

const b = new Float64Array([
  2, 4, 8, 6
])

計算の為に、適当な値いれておいた。

初期近似解が真に近い値で初期化されるほど探索が減り高速に計算が終わることになる。

function CGSolver (A, x, b) {
  const n = x.length

  let r = new Float64Array(n) // 残差
  let p = new Float64Array(n) // 基底ベクトル

  // ↓ 初期近似解による残差
  for (let i = 0; i < n; ++i) {
    let Ax = 0.0
    for (let j = 0; j < n; ++j) {
      Ax += A[j * n + i] * x[j]
    }
    r[i] = b[i] - Ax
    p[i] = r[i]
  }

  let rr0 = dot(r, r)

  let Ap = new Float64Array(n)

  const limit = 40

  for (let k = 0; k < limit; ++k) {

    // Ap = (A, p)
    for (let i = 0; i < n; ++i) {
      let res = 0
      for (let j = 0; j < n; ++j) {
        res += A[i * n + j] * p[j]
      }
      Ap[i] = res
    }

    // ↓ alpha = (r, r) / (p, AP)
    const alpha = rr0 / dot(p, Ap)

    // ↓ 近似解と残差の修正
    for (let i = 0; i < n; ++i) {
      x[i] += alpha * p[i]
      r[i] -= alpha * Ap[i]
    }

    const rr1 = dot(r, r)

    // ↓ ベクトルノルム Σ|x[i]| < eps
    let norm = 0
    for (let i = 0; i < n; ++i) {
      norm += Math.abs(r[i])
    }
    if (norm < 1.0e-6) break

    // ↓ beta = (r, r) / (r', r')
    const beta = rr1 / rr0
    for (let i = 0; i < n; ++i) {
      p[i] = r[i] + beta * p[i]
    }

    rr0 = rr1
  }
}

function dot (x, y) {
  let res = 0
  for (let i = 0; i < x.length; ++i) {
    res += x[i] * y[i]
  }
  return res
}

反復の最大数は40に設定、収束と判定された時点で計算は終了。

ガウスの消去法を用いて誤差を確認する。

Node
const A = new Float64Array([
  1, 2, 4, 6,
  2, 3, 2, 8,
  4, 2, 2, 4,
  6, 8, 4, 1
])

const x = new Float64Array([
  0, 0, 0, 0
])

const b = new Float64Array([
  2, 4, 8, 6
])

const _A = A.slice()
const _x = x.slice()
const _b = b.slice()

CGMethod(A, x, b)

gaussianElimination(_A, _x, _b)

for (let i = 0, len = x.length; i < len; ++i) {
  const d = Math.abs(x[i] - _x[i])
  console.log(d)
}

function gaussianElimination (A, x, b) {
  const n = x.length

  let c = 0
  for (let i = 0; i < n - 1; i++) {
    for (let j = i + 1; j < n; j++) {
      c = A[j * n + i] / A[i * n + i]
      for (let k = i + 1; k < n; k++) {
        A[j * n + k] -= c * A[i * n + k]
      }
      b[j] -= c * b[i]
    }
  }
  for (let i = n - 1; i >= 0; i--) {
    let sum = 0.0
    for (let j = i + 1; j < n; j++) {
      sum += A[i * n + j] * x[j]
    }
    x[i] = (b[i] - sum) / A[i * n + i]
  }
}

計算の誤差は以下のようになる。

|x[0] - _x[0]| = 2.6645352591003757e-15
|x[1] - _x[1]| = 4.3298697960381105e-15
|x[2] - _x[2]| = 1.0894063429134349e-15
|x[3] - _x[3]| = 4.2188474935755950e-15

3回の計算速度は以下のようになる。

ICCGSolver 0.770ms 0.673ms 0.670ms
gaussianElimination 0.284ms 0.250ms 0.255ms

遅くなっているのは、行列の次元が低いことや初期近似解が0であることなど。

不完全コレスキー分解

コレスキー分解(コレスキーぶんかい)とは、正定値エルミート行列Aを下三角行列LとLの共役転置L*との積に分解すること。
= ja.wikipedia.org

Node
function incompleteCholeskyFactorization (A, L, d) {
  const n = d.length

  L[0] = A[0]
  d[0] = 1.0 / L[0]
  for (let i = 1; i < n; ++i) {
    for (let j = 0; j <= i; ++j) {
      if (Math.abs(A[j * n + i]) < 1.0e-10) continue
      let lld = A[j * n + i]
      for (let k = 0; k < j; ++k) {
        lld -= L[k * n + i] * L[k * n + j] * d[k]
      }
      L[j * n + i] = lld
    }
    d[i] = 1.0 / L[i * n + i]
  }
}

分解された対角成分と下三角行列を生成するのみ。

Node
const A = new Float64Array([
  1, 2, 4, 6,
  2, 3, 2, 8,
  4, 2, 2, 4,
  6, 8, 4, 1
])
let d = new Float64Array(n) // 対角成分
let L = new Float64Array(n * n) // 下三角行列

incompleteCholeskyFactorization(A, L, d)

不完全コレスキー分解付共役勾配法

(LDL^t)^-1rの内積を求める処理は関数ICRes。

function ICCGSolver (A, x, b) {
  const n = x.length

  let r = new Float64Array(n) // 残差
  let p = new Float64Array(n) // 基底ベクトル

  let r2 = new Float64Array(n)

  let d = new Float64Array(n) // 対角成分
  let L = new Float64Array(n * n) // 下三角行列

  incompleteCholeskyFactorization(A, L, d)

  // ↓ 初期近似解による残差
  for (let i = 0; i < n; ++i) {
    let Ax = 0.0
    for (let j = 0; j < n; ++j) {
      Ax += A[j * n + i] * x[j]
    }
    r[i] = b[i] - Ax
  }

  // ↓ ((LDL^t)^-1, r)
  ICRes(L, d, r, p, n)

  let rr0 = dot(r, p, n)

  let Ap = new Float64Array(n)

  const limit = 40

  for (let k = 0; k < limit; ++k) {

    // Ap = (A, p)
    for (let i = 0; i < n; ++i) {
      let res = 0
      for (let j = 0; j < n; ++j) {
        res += A[i * n + j] * p[j]
      }
      Ap[i] = res
    }

    // ↓ alpha = (r,r) / (p,Ap)
    let alpha = rr0 / dot(p, Ap, n)

    // ↓ 近似解と残差の修正
    for (let i = 0; i < n; ++i) {
      x[i] += alpha * p[i]
      r[i] -= alpha * Ap[i]
    }

    ICRes(L, d, r, r2)

    let rr1 = dot(r, r2, n)

    // ↓ ベクトルノルム Σ|x[i]| < eps
    let norm = 0
    for (let i = 0; i < n; ++i) {
      norm += Math.abs(r[i])
    }
    if (norm < 1.0e-6) break

    // ↓ beta = (r, r) / (r', r')
    let beta = rr1 / rr0
    for (let i = 0; i < n; ++i) {
      p[i] = r2[i] + beta * p[i]
    }

    rr0 = rr1
  }
}

function dot (x, y) {
  let res = 0
  for (let i = 0; i < x.length; ++i) {
    res += x[i] * y[i]
  }
  return res
}

function incompleteCholeskyFactorization (A, L, d) {
  const n = d.length

  L[0] = A[0]
  d[0] = 1.0 / L[0]
  for (let i = 1; i < n; ++i) {
    for (let j = 0; j <= i; ++j) {
      if (Math.abs(A[j * n + i]) < 1.0e-10) continue
      let lld = A[j * n + i]
      for (let k = 0; k < j; ++k) {
        lld -= L[k * n + i] * L[k * n + j] * d[k]
      }
      L[j * n + i] = lld
    }
    d[i] = 1.0 / L[i * n + i]
  }
}

function ICRes (L, d, r, u, n) {
  let y = new Float64Array(n)

  for (let i = 0; i < n; ++i) {
    let rly = r[i]
    for (let j = 0; j < i; ++j) {
      rly -= L[j * n + i] * y[j]
    }
    y[i] = rly / L[i * n + i]
  }

  for (let i = n - 1; i >= 0; --i) {
    let lu = 0.0
    for (let j = i + 1; j < n; ++j) {
      lu += L[i * n + j] * u[j]
    }
    u[i] = y[i] - d[i] * lu
  }
}

計算の誤差は以下のようになる。

|x[0] - _x[0]| = 8.8817841970012520e-16
|x[1] - _x[1]| = 4.4408920985006260e-16
|x[2] - _x[2]| = 2.7755575615628914e-17
|x[3] - _x[3]| = 1.1102230246251565e-16

3回の計算速度は以下のようになる。

ICCGSolver 1.040ms 1.191ms 1.384ms
gaussianElimination 0.251ms 0.262ms 0.286ms

不完全コレスキー分解で高速化される行列は限定されている。