VBAで行列計算を行う様々なコードを実装しました。足し算、引き算、掛け算、転置、単位行列作成、掃き出し法、逆行列、水平結合、疑似逆行列、行列式、型変換、セルからのデータ読み込みと書き出しです。個人的には、もうVBAを触りたくないので封印したい。笑

matrix.vba

Function matrix_plus(m1, m2)
  ' 行列同士の足し算
  ans = m1 ' 同じ大きさのオブジェクトを生成
  For i = 0 To UBound(m1)
    For j = 0 To UBound(m1(0))
      ans(i)(j) = m1(i)(j) + m2(i)(j)
    Next j
  Next i
  matrix_plus = ans
End Function




Function matrix_minus(m1, m2)
  ' 行列同士の引き算
  ans = m1 ' 同じ大きさのオブジェクトを生成
  For i = 0 To UBound(m1)
    For j = 0 To UBound(m1(0))
      ans(i)(j) = m1(i)(j) - m2(i)(j)
    Next j
  Next i
  matrix_minus = ans
End Function


Function create_matrix(row_size, col_size)
    ' 任意サイズの行列を作成する
    Dim ans, row As Variant
    
    ans = Array()
    ReDim ans(row_size - 1)
    For i = 0 To row_size - 1
        row = Array() ' 新しいオブジェクトのインスタンスが代入される
        ReDim row(col_size - 1)
        ans(i) = row
    Next i
    
    create_matrix = ans
End Function


Function matrix_cross(m1, m2)
  ' 行列同士の掛け算
  ans = create_matrix(UBound(m1) + 1, UBound(m2(0)) + 1)
  For i = 0 To UBound(ans)
    For j = 0 To UBound(ans(0))
      sum_ = 0
      For k = 0 To UBound(m1(0))
        sum_ = sum_ + m1(i)(k) * m2(k)(j)
      Next k
      ans(i)(j) = sum_
    Next j
  Next i
  matrix_cross = ans
End Function


Function matrix_t(m)
  ' 行列の転置
  ans = create_matrix(UBound(m(0)) + 1, UBound(m) + 1)
  For i = 0 To UBound(ans)
    For j = 0 To UBound(ans(0))
      ans(i)(j) = m(j)(i)
    Next j
  Next i
  matrix_t = ans
End Function


Function create_identity_matrix(size)
  ' 単位行列を作成する
  ans = create_matrix(size, size)
 
  For i = 0 To size - 1
    For j = 0 To size - 1
      If i = j Then
        ans(i)(j) = 1
      Else
        ans(i)(j) = 0
      End If
    Next j
  Next i
  create_identity_matrix = ans
End Function



'==========================================================================
' ベクトル（1次元配列やArray）に関する変数の型変換
Function vector2matrix(vect)
  ' ベクトルを行列に変換する（列ベクトル扱いにする）
  matrix = Array(vect) ' このままでは行ベクトル扱い
  vector2matrix = matrix_t(matrix) ' これで行列の中で列ベクトル扱いになる
End Function


Sub matrix_test()
    row1 = Array(11, 12, 13)
    row2 = Array(21, 22, 23)
    row3 = Array(31, 32, 33)
    temp_matrix1 = Array(row1, row2, row3)
    
    result = matrix_plus(temp_matrix1, temp_matrix1)
    
    foo = create_matrix(3, 4)
    foo(0)(0) = 100
    
    row1 = Array(11, 12, 13, 14)
    row2 = Array(21, 22, 23, 24)
    row3 = Array(31, 32, 33, 34)
    temp_matrix2 = Array(row1, row2, row3)
    
    row1 = Array(11, 12, 13)
    row2 = Array(21, 22, 23)
    row3 = Array(31, 32, 33)
    row4 = Array(41, 42, 43)
    temp_matrix3 = Array(row1, row2, row3, row4)
    
    result = matrix_cross(temp_matrix2, temp_matrix3)
    result = matrix_t(result)
    
    result = create_identity_matrix(10)
    
    ' ベクトルと行列の積
    vect = Array(1, 2, 3)
    m = vector2matrix(vect)
    result = matrix_cross(temp_matrix1, m)
    result = matrix_cross(m, temp_matrix1) ' エラーになってほしいがならない。
End Sub



'==========================================================================
' 掃き出し法関連
Function search_max_coef(search_col, matrix_data)
 ' search_col列の最大値を持つ行を返す。ただし、search_col行以降を走査する。
  row_at_max = -1
  max_val = 0
  For row = search_col To UBound(matrix_data)
    Value = Abs(matrix_data(row)(search_col))
    If Value > max_val Then
      max_val = Value
      row_at_max = row
    End If
  Next row
  search_max_coef = row_at_max
End Function

Function switch_row(row_index1, row_index2, matrix_data)
 ' row_index1行目とrow_index2行目を入れ替える
  row1 = matrix_data(row_index1)
  row2 = matrix_data(row_index2)
  matrix_data(row_index1) = row2
  matrix_data(row_index2) = row1
  switch_row = matrix_data
End Function

Function matrix_append(matrix_A_, vector_Y_)
    ' 行列の最終列にvector_Y_の要素をコピーする。実用には要素数のチェックを入れたほうが良い。
    new_matrix = matrix_A_
    For i = 0 To UBound(new_matrix)
        Dim row_data As Variant
        row_data = new_matrix(i)
        ReDim Preserve row_data(UBound(row_data) + 1) ' 内容保持したまま拡張
        row_data(UBound(row_data)) = vector_Y_(i)
        new_matrix(i) = row_data
    Next i
    matrix_append = new_matrix
End Function

Function get_x_0(vector_Y_, matrix_A_)
  ' 掃き出し法により解を求める
  vector_Y = vector_Y_ ' Variant型ならDeep copyされる
  matrix_a = matrix_A_
  For row = 0 To UBound(vector_Y)
    ' row行の係数を(row行, row列)の係数で割る
    coef = matrix_a(row)(row)
    vector_Y(row) = vector_Y(row) / coef
    For col = row To UBound(vector_Y)
      matrix_a(row)(col) = matrix_a(row)(col) / coef
    Next col
    ' row行に(row2行, row列)の係数を掛けたものをrow2行から引く。ただしrow=row2は除く。
    For row2 = 0 To UBound(vector_Y)
      If row <> row2 Then
        coef = matrix_a(row2)(row)
        vector_Y(row2) = vector_Y(row2) - vector_Y(row) * coef
        For col = row To UBound(vector_Y)
          matrix_a(row2)(col) = matrix_a(row2)(col) - matrix_a(row)(col) * coef
        Next col
      End If
    Next row2
  Next row
  get_x_0 = vector_Y
End Function


Function get_x_1(vector_Y_, matrix_A_)
  ' 掃き出し法により解を求める
  matrix_data = matrix_append(matrix_A_, vector_Y_)
  For row = 0 To UBound(matrix_data)
    ' row行の係数を(row行, row列)の係数で割る
    coef = matrix_data(row)(row)
    For col = row To UBound(matrix_data(0))
      matrix_data(row)(col) = matrix_data(row)(col) / coef
    Next col

    ' row行に(row2行, row列)の係数を掛けたものをrow2行から引く。ただしrow=row2は除く。
    For row2 = 0 To UBound(matrix_data)
      If row <> row2 Then
        coef = matrix_data(row2)(row)
        For col = row To UBound(matrix_data(0))
          matrix_data(row2)(col) = matrix_data(row2)(col) - matrix_data(row)(col) * coef
        Next col
      End If
    Next row2
  Next row

  ' 解を取り出す(なんて面倒なんだ）
  Dim ans As Variant
  ans = Array()
  ReDim ans(UBound(matrix_data))
  For i = 0 To UBound(matrix_data)
    ans(i) = matrix_data(i)(UBound(matrix_data(0)))
  Next i
  
  get_x_1 = ans
End Function

Function get_x_2(vector_Y_, matrix_A_)
  ' 掃き出し法により解を求める
  matrix_data = matrix_append(matrix_A_, vector_Y_)
  For row = 0 To UBound(matrix_data)
    ' ピボット操作
    switch_row_index = search_max_coef(row, matrix_data)
    matrix_data = switch_row(row, switch_row_index, matrix_data)

    ' row行の係数を(row行, row列)の係数で割る
    coef = matrix_data(row)(row)
    For col = row To UBound(matrix_data(0))
      matrix_data(row)(col) = matrix_data(row)(col) / coef
    Next col

    ' row行に(row2行, row列)の係数を掛けたものをrow2行から引く。ただしrow=row2は除く。
    For row2 = 0 To UBound(matrix_data)
      If row <> row2 Then
        coef = matrix_data(row2)(row)
        For col = row To UBound(matrix_data(0))
          matrix_data(row2)(col) = matrix_data(row2)(col) - matrix_data(row)(col) * coef
        Next col
      End If
    Next row2
  Next row

  ' 解を取り出す(なんて面倒なんだ）
  Dim ans As Variant
  ans = Array()
  ReDim ans(UBound(matrix_data))
  For i = 0 To UBound(matrix_data)
    ans(i) = matrix_data(i)(UBound(matrix_data(0)))
  Next i
  
  get_x_2 = ans
End Function


Function matrix_join(matrix_A_, matrix_B_)
    ' 行列の右にmatrix_B_の要素をコピーする。実用には要素数のチェックを入れたほうが良い。
    ' この関数単体ではあまり役に立たないが、掃き出し法の関数を整理すると逆行列を求める際に活きてくる。
    new_matrix = matrix_A_ ' deep copy. 行数はOKだがまだ列数は半分。
    For i = 0 To UBound(new_matrix)
        Dim row_data As Variant
        row_data = new_matrix(i)
        w = UBound(row_data) + 1
        ReDim Preserve row_data(w + UBound(matrix_B_(0)))  ' 内容保持したまま拡張
        For j = 0 To UBound(matrix_B_(0))
            row_data(w + j) = matrix_B_(i)(j)
        Next j
        new_matrix(i) = row_data
    Next i
    matrix_join = new_matrix
End Function


Sub hakidashi_test()
    row1 = Array(1, -2, -2, 2)
    row2 = Array(2, -2, -3, 3)
    row3 = Array(-1, 6, 3, -2)
    row4 = Array(1, 4, 0, -1)
    temp_matrix = Array(row1, row2, row3, row4)
    y = Array(5, 10, 2, -10)
    result = get_x_0(y, temp_matrix)
    result = get_x_1(y, temp_matrix)
    result = get_x_2(y, temp_matrix)
        
    temp_matrix = Array(row1, row2, row3, row4)
    Cells(5, 1) = search_max_coef(0, temp_matrix)
    Cells(6, 1) = search_max_coef(1, temp_matrix)
    Cells(7, 1) = search_max_coef(2, temp_matrix)
    Cells(8, 1) = search_max_coef(3, temp_matrix)
    
    temp_matrix = Array(row1, row2, row3, row4)
    temp_matrix2 = switch_row(0, 1, temp_matrix)
    
    result = matrix_join(temp_matrix, temp_matrix)
End Sub


'==========================================================================
Function matrix_inverse(matrix)
    ' 逆行列を求める
    Dim delta As Variant
    
    delta = Array()
    ReDim delta(UBound(matrix))
    inverse = matrix ' コピーはしているが、ただ同じサイズの行列が欲しいだけ
    For i = 0 To UBound(matrix)
        For k = 0 To UBound(delta)
            If k = i Then
                delta(k) = 1 ' 残りは0で初期化されることを期待
            Else
                delta(k) = 0
            End If
        Next k
        b = get_x_2(delta, matrix)
        inverse(i) = b
    Next i
    inverse = matrix_t(inverse)
    
    matrix_inverse = inverse
End Function

Sub inverse_test()
    row1 = Array(2, 5)
    row2 = Array(1, 3)
    temp_matrix = Array(row1, row2)
    
    result = matrix_inverse(temp_matrix)
    result = matrix_cross(temp_matrix, result) ' 単位行列になれば成功
End Sub




'==========================================================================
Function matrix_pseudo_inverse(m)
    ' 疑似逆行列を返す
    matrix_pseudo_inverse = matrix_cross(matrix_inverse(matrix_cross(matrix_t(m), m)), matrix_t(m))
End Function

Sub pseudo_inverse_test()
    row1 = Array(2, 5)
    row2 = Array(1, 3)
    row3 = Array(4, 7)
    row4 = Array(8, 2)
    temp_matrix = Array(row1, row2, row3, row4)
    
    result = matrix_pseudo_inverse(temp_matrix)
    result = matrix_cross(result, temp_matrix) ' 単位行列になれば成功。掛ける順番（引数の順番）が重要。
End Sub




'==========================================================================
' 行列式を求める
' 参考:http://thira.plavox.info/blog/2008/06/_c.html
' 上の参考Webサイトのプログラムは行列式が0の場合に計算できないので、少し改造する。
Function matrix_det(m_)
    ' 引数が正方行列であることをまずチェックした方が良いが、ここでは省略する
    m = m_
    n = UBound(m) '配列の次数-1
    
    For i = 0 To n
        ' i列が全部0にだったら行列式は0
        all_zero = True
        For k = 0 To n
            If Abs(m(k)(i)) > 0.0000000001 Then
                all_zero = False
            End If
        Next k
        If all_zero = True Then ' 関数から脱出
            matrix_det = 0
            Exit Function
        End If
                
                
        ' ピボット操作。i行以降でi列が最大の行とi行を入れ替える
        switch_row_index = search_max_coef(i, m)
        m = switch_row(i, switch_row_index, m)
        
        For j = 0 To n
            If i < j Then
                ' 上三角行列を作る
                buf = m(j)(i) / m(i)(i)
                For k = 0 To n
                    m(j)(k) = m(j)(k) - m(i)(k) * buf
                Next k
                
                ' j行が全部0になったら行列式は0
                all_zero = True
                For k = 0 To n
                    If Abs(m(j)(k)) > 0.0000000001 Then
                        all_zero = False
                    End If
                Next k
                If all_zero = True Then ' 関数から脱出
                    matrix_det = 0
                    Exit Function
                End If
            End If
        Next j
    Next i
    
    det = 1#
    For i = 0 To n
        det = det * m(i)(i)
    Next i
    
    matrix_det = det
End Function


Sub det_test()
    a1 = Array(2, -2, 4, 2)
    a2 = Array(2, -1, 6, 3)
    a3 = Array(3, -2, 12, 12)
    a4 = Array(-1, 3, -4, 4)
    a = Array(a1, a2, a3, a4)
    det = matrix_det(a) '=>120
    
    a1 = Array(10, 11, 12, 13) ' 連立方程式の係数行列。連立方程式は数学的には解なしとなる問題
    a2 = Array(13, 15, 14, 13)
    a3 = Array(10, 11, 12, 13)
    a4 = Array(14, 15, 16, 17)
    a = Array(a1, a2, a3, a4)
    det = matrix_det(a) '=>0
    
    a1 = Array(3, 2, -3, 1) ' 連立方程式の係数行列。連立方程式は数学的には無数の解が得られる問題
    a2 = Array(4, 3, -4, 2)
    a3 = Array(2, 1, 2, -1)
    a4 = Array(6, 3, 2, -2)
    a = Array(a1, a2, a3, a4)
    det = matrix_det(a) '=>0
    
    a1 = Array(2, 1, 3) ' 連立方程式の係数行列。連立方程式は数学的には解がある問題
    a2 = Array(3, -1, 4)
    a3 = Array(1, 3, -2)
    a = Array(a1, a2, a3)
    det = matrix_det(a) '=>20
End Sub




'==========================================================================
' 行列に関する変数の型変換
Function change2DtoArray(matrix2d)
    ' 普通の2次元配列をArrayの入れ子構造に変換する
    m = UBound(matrix2d, 1) + 1 ' 行数
    n = UBound(matrix2d, 2) + 1 ' 列数

    matrix = create_matrix(m, n)

    For i = 0 To m - 1
        For j = 0 To n - 1
            matrix(i)(j) = matrix2d(i, j)
        Next j
    Next i
    change2DtoArray = matrix
End Function


Function changeArrayTo2D(matrix_array)
    ' Arrayの入れ子構造で表した行列を普通の2次元配列での表現に変換する
    m = UBound(matrix_array)    ' 行数-1
    n = UBound(matrix_array(0)) ' 列数-1

    Dim matrix2d() As Double
    ReDim matrix2d(m, n)

    For i = 0 To m
        For j = 0 To n
            matrix2d(i, j) = matrix_array(i)(j)
        Next j
    Next i
    changeArrayTo2D = matrix2d
End Function

Sub type_change_test()
    Dim a(2, 2) As Double
    a(0, 0) = 11
    a(0, 1) = 12
    a(0, 2) = 13
    a(1, 0) = 21
    a(1, 1) = 22
    a(1, 2) = 23
    a(2, 0) = 31
    a(2, 1) = 32
    a(2, 2) = 33
    
    Dim b(2) As Double
    b(0) = 1
    b(1) = 2
    b(2) = 3
    
    Dim c(2) As Double
    c(0) = 1
    c(1) = 2
    c(2) = 3
    
    a_ = change2DtoArray(a)
    a__ = changeArrayTo2D(a_) ' aと同一であればOK
    b_ = vector2matrix(b)
    c_ = vector2matrix(c)
End Sub



'==========================================================================
' data read from sheet and write to sheet
Function read_matrix_from_sheet(row_origin, col_origin, row_size, col_size)
    ' シートから原点と縦横のサイズを指定し、2次元状にデータを取得する
    ' Range関数の代わりに使える。
    ' col_origin is origin of col (upper left of matrix on sheet)
    ' row_origin is origin of row (upper left of matrix on sheet)
    ' col_size is matrix col size
    ' row_size is matrix row size

    Dim a() As Variant
    ReDim a(row_size - 1, col_size - 1)
    
    For i = 0 To UBound(a, 1)
        For j = 0 To UBound(a, 2)
            a(i, j) = Cells(row_origin + i, col_origin + j)
        Next j
    Next i
    read_matrix_from_sheet = a

End Function

Function change_zero_start(m)
    ' インデックスが1から始まる2次元配列を0から始まる2次元配列に変換する
    ' Range関数で取得した2次元配列を変換することを想定している。
    ' m is a 2d-array as matrix. its index starts from 1.
    ' this function returns a 0-start-matrix.

    Dim a() As Variant
    ReDim a(UBound(m, 1) - 1, UBound(m, 2) - 1)
    
    For i = LBound(m, 1) To UBound(m, 1)
        For j = LBound(m, 2) To UBound(m, 2)
            a(i - 1, j - 1) = m(i, j)
        Next j
    Next i
    change_zero_start = a
End Function

Sub read_write_test()
    ' read from cells
    ' セルから読み込んだり書き込んだりするサンプル
    matrix_a = read_matrix_from_sheet(1, 1, 3, 4)
    
    hoge = Range("b2:e4") ' index is 1-start.
    fuga = change_zero_start(hoge)


    ' write to cells
    row1 = Array(11, 12, 13) ' 行ベクトルのつもり
    row2 = Array(21, 22, 23)
    row3 = Array(31, 32, 33)
    m = Array(row1, row2, row3)
    'Range("A1:F30") = m  ' この書き方は型の問題でうまくいかない
    Range("A1:F30") = changeArrayTo2D(m) ' 普通の2次元配列なら、一度に書き出せるようだ。
End Sub

9行目がmatrix_plusでなくmatrix_minusになっています。

ありがとうございます。修正しました。