Numerical Recipes in CのCholesky分解ルーチンをExcel VBAに移植

ニューメリカルレシピをExcel　VBAに移植するシリーズ。今回はコレスキー分解。

これは非常に簡単に、

Sub choldc(a() As Double, n As Long, p() As Double)
     Dim i As Long, j As Long, k As Long
     Dim sum As Double
    
     For i = 1 To n
        For j = i To n
            sum = a(i, j)
            For k = i - 1 To 1 Step -1
                sum = sum - a(i, k) * a(j, k)
            Next k
            If (i = j) Then
                If (sum <= 0#) Then
                    MsgBox "choldc failed"
                End If
                p(i) = Sqr(sum)
            Else
                a(j, i) = sum / p(i)
            End If
        Next j
    Next i

End Sub

と書ける。

使い方は、aを2次元配列で行列を表し、正の値を持つ対称行列とする。で上三角行列だけを入力しておけばOKで、下三角行列にコレスキー分解した値が入る。ただし、対角成分は一次元配列のpに入る。

例えば5x5のHilbert行列 ( aij = 1/(i+j-1) )をコレスキー分解するには、

Option Explicit

Private Sub CommandButton1_Click()
    Dim n As Long
    Dim a() As Double, p() As Double, b() As Double
    Dim i As Long, j As Long
   
    n = 5
    ReDim a(n, n)
    ReDim b(n, n)
    ReDim p(n)
   

   
    For i = 1 To n
        For j = i To n
            a(i, j) = 1# / (i + j - 1)
            Worksheets("sheet1").Cells(i + 1, j + 1) = a(i, j)
        Next j
    Next i

    Call choldc(a, n, p)
   
   
   
    For j = 1 To n
        For i = j To n
            Worksheets("sheet1").Cells(i + 1, j + 7) = a(i, j)
        Next i
    Next j
    For i = 1 To n

            Worksheets("sheet1").Cells(i + 1, n + 9) = p(i)

    Next i
End Sub

で呼び出す。Maximaで計算した結果がこちら。

でばっちり合ってます。

さて、昨日メルセンヌツイスタで正規乱数を計算できるようになって、本日コレスキー分解ができるようになったということは？（続く）。

----

Excel　VBAにNumerical Recipesのルーチンを移植するシリーズ：

三重対角行列

https://sci.tea-nifty.com/blog/2012/07/numerical-recip.html

LU分解

https://sci.tea-nifty.com/blog/2010/07/excel-vbanumeri.html