行列の掛け算をライブラリを使わず自力で実装する

Question

事情があってライブラリを使わず最初からPythonに入っているもので行列の掛け算を実装しなければいけません。

一応できたのはできたのですが、コードの正当性の証明などどうやればいいかわからないです。
行列の掛け算をやってくれるサイトがあるのでそれの答えと照らし合わせていますが正しいか不安です。
そもそもあっているのかわからないです。
よろしくお願いします。
参考

#行列UとVを掛け算する
def matrix_multiplication(U,V):
    answer=[]
    u_colmun_length=len(U[0])
    v_row_lengthh=len(V)
    v_column_length=len(V[0])
    #Uの列数とVの行数が同じじゃないなら計算できない
    if u_colmun_length != v_row_lengthh:
        print("Could not compute")
        return
    else:
        #Uの行列から１行を取り出す
        for u in U:
            temp=[]        
            colmun_index=0
            #Vの行列から一列を取り出す
            for _ in range(v_column_length):
                a = [x[colmun_index] for x in V]
                el_index=0         
                temp_multiple=0
                #Uの行とVの列を掛け算し足し合わせる
                for _ in range(u_colmun_length):   
                    temp_multiple+=u[el_index]*a[el_index] 
                    el_index+=1         
                #行列の行を作成
                temp.append(temp_multiple)                
                colmun_index+=1   
            #行を追加
            answer.append(temp)
    print(answer)

matrix_multiplication([[2,3],[3,5]],[[2,1]]) 
matrix_multiplication([[2,3]],[[2],[3]])
matrix_multiplication([[2,3],[2,10]],[[3,5],[3,3]]) 
matrix_multiplication([[2,3,3],[2,10,3]],[[3,5,2],[3,3,1],[2,4,5]])       

"""
Result

Could not compute
[[13]]
[[15, 19], [36, 40]]
[[21, 31, 22], [42, 52, 29]]
"""

Answer 1

コードの正当性の証明などどうやればいいかわからないです

普通は、テストケースを作り、テストに全て（すべてでない場合もある）合格すれば良しとします。
このような計算の場合、事前に答えがわかっている値を何組か用意して、確認します。

というので回答になりますか？他の何かを求めていますか？

Answer 2

10000 回ランダムに生成した行列で試しましたが、すべて numpy の答えと一致しましたので、実装は合っていると思います。

# 行列UとVを掛け算する
def matrix_multiplication(U, V):
    answer = []
    u_colmun_length = len(U[0])
    v_row_lengthh = len(V)
    v_column_length = len(V[0])
    # Uの列数とVの行数が同じじゃないなら計算できない
    if u_colmun_length != v_row_lengthh:
        print("Could not compute")
        return None
    else:
        # Uの行列から１行を取り出す
        for u in U:
            temp = []
            colmun_index = 0
            # Vの行列から一列を取り出す
            for _ in range(v_column_length):
                a = [x[colmun_index] for x in V]
                el_index = 0
                temp_multiple = 0
                # Uの行とVの列を掛け算し足し合わせる
                for _ in range(u_colmun_length):
                    temp_multiple += u[el_index] * a[el_index]
                    el_index += 1
                # 行列の行を作成
                temp.append(temp_multiple)
                colmun_index += 1
            # 行を追加
            answer.append(temp)
        return answer
########################################################
# ここまで質問者さんのコード
########################################################
    
from pprint import pprint
import numpy as np

# 計算できるパターンで答えがあっているかどうか
np.random.seed(42)
num_tests = 10000
num_ok = 0

for i in range(num_tests):
    # 行列 A, B の列数、行数をランダムに生成する。
    N, M, L = np.random.randint(1, 10, 3)
    # 行列 A, B の値をランダムに生成する。
    A = np.random.randint(1, 10, size=(N, M))
    B = np.random.randint(1, 10, size=(M, L))

    # numpy の計算結果
    ret1 = A @ B

    # 自作した関数の計算結果
    ret2 = matrix_multiplication(A, B)

    if not np.array_equal(ret1, ret2):
        print("=== incorrect case found ===")
        print(f"A {A.shape}\n{A}")
        print(f"B {B.shape}\n{B}")
        print(f"numpy {ret1.shape}:\n{ret1}")
        print(f"matrix_multiplication:\n{ret2}")
    else:
        num_ok += 1
        
print(f"{num_ok} / {num_tests}")  # 10000 / 10000

補足

A: (N, M) の行列, B: (M, L) の行列としたとき、
行列積 C は (N, L) の行列でその (i, j) 成分は

C_{ij} = \sum_k^M A_{ik} B_{kj}

なので、定義に従い書くと以下のようになる。

def matrix_multiplication(A, B):
    A_rows, A_cols = len(A), len(A[0])  # A の行数及び列数
    B_rows, B_cols = len(B), len(B[0])  # B の行数及び列数

    if A_cols != B_rows:
        print("Could not compute")
        return None

    # (A_rows, B_cols) の0で初期化した行列を作成する。
    ret = [[0 for j in range(B_cols)] for i in range(A_rows)]
    
    for i in range(A_rows):
        for j in range(B_cols):
            for k in range(A_cols):
                ret[i][j] += A[i][k] * B[k][j]

    return ret

追記

↓これは[[0]]と出るのですが二次元配列ですか。
ret = [[0 for j in range(B_cols)] for i in range(A_rows)]
ここにインデックスをあらかじめ指定して代入するということですかね？
[[0]]と表示されているのにIndex out of rangeみたいなエラーが出ないので不思議でした

例えば、(2, 2) と (2, 3) の行列積であれば、結果は (2, 3) の行列となるので、(2, 3) の値がすべて0の行列を Python の内包記法で作成しています。

pythonの内包表記を少し詳しく - Qiita

ret = [[0 for j in range(3)] for i in range(2)]
print(ret)
# [[0, 0, 0],
#  [0, 0, 0]]

具体的な例で各出力がどうなるかをコメントに記載したコード

A = [[1, 2],
     [3, 4]]
B = [[5, 6, 7],
     [8, 9, 10]]

A_rows, A_cols = len(A), len(A[0])  # A の行数及び列数
print(f"({A_rows}, {A_cols})")  # (2, 2)
B_rows, B_cols = len(B), len(B[0])  # B の行数及び列数
print(f"({B_rows}, {B_cols})")  # (3, 2)

if A_cols != B_rows:
    print("Could not compute")
    exit(1)

# 0で初期化した 2x3 行列を作成する。
ret = [[0 for j in range(B_cols)] for i in range(A_rows)]
print(ret)
# [[0, 0, 0],
#  [0, 0, 0]]

for i in range(A_rows):
    for j in range(B_cols):
        for k in range(A_cols):
            ret[i][j] += A[i][k] * B[k][j]

あと↓ですがかけたものを足しているということでいいですか？勉強不足のせいか読み方が分かりません・・・。
C_{ij} = \sum_k^M A_{ik} B_{kj}

行列積の定義です。Tex 記法が teratail だとコンパイルされないので、tex のままになっていますが、以下の行列積の定義を書きました。

Comments

[eggpol]

回答ありがとうございます。
やはり確認するにはテストケースを作って確かめるのですね。
コードを拝見いたしましたがこんなにスマートに行列の掛け算ができるのですね。

↓これは[[0]]と出るのですが二次元配列ですか。
ret = [[0 for j in range(B_cols)] for i in range(A_rows)]
ここにインデックスをあらかじめ指定して代入するということですかね？
[[0]]と表示されているのにIndex out of rangeみたいなエラーが出ないので不思議でした

あと↓ですがかけたものを足しているということでいいですか？勉強不足のせいか読み方が分かりません・・・。
C_{ij} = \sum_k^M A_{ik} B_{kj}

[tiitoi]

回答に追記しました。
まだ不明点がございましたら、コメントしてください。

Answer 3

正しいかどうかを確認するときにライブラリを使ったらいいではないですか？