cms.matrix.matrix のソースコード

"""
概要: NumPyとSciPyを用いた基本的な行列演算の例を示すスクリプト。

詳細説明:
このスクリプトは、PythonにおけるNumPyおよびSciPyライブラリを使用して、様々な行列演算を実行する方法を実演します。
行列の初期化、内積、外積、逆行列、行列式、固有値問題、連立一次方程式の解法、LU分解、
コレスキー分解、QR分解などの具体的なコード例を提示します。
各演算の結果は標準出力に表示され、線形代数の基本的な操作を理解するのに役立ちます。

関連リンク: :doc:`matrix_usage`
"""
import sys
import numpy as np
import scipy as sp
import numpy.linalg as npalg
import scipy.linalg as spalg

# Better not use matrix type, instead use ndarray
# Q: real orthogonal matrix
# R: upper-triangle matrix
# S: posiive definite symmetric matrix

#===================
# Global parameters (Initial Data)
#===================
A  = np.array([[2, 5], 
               [1, 3]])
B  = np.array([[-1], 
               [ 1]])
S  = np.array([[ 2.0, -1.0]
             , [-1.0,  5.0]])
V1 = np.array([1, 1, 2])
V2 = np.array([2, 1, 3])



def main():
    """
    概要: 行列演算の様々な機能を実行し、その結果を標準出力に表示する。

    詳細説明:
    この関数は、NumPyおよびSciPyの線形代数機能を用いて、事前に定義された行列 A, B, S とベクトル V1, V2
    に対して、内積、外積、逆行列、行列式、固有値計算、連立方程式の解法、LU分解、コレスキー分解、
    QR分解などの基本的な行列演算を行います。
    各演算の結果はコンソールに出力され、線形代数の概念をPythonでどのように実装するかを示します。
    """
    print("")
    print("A:\n", A)
    print("B:\n", B)
    print("S:\n", S)
    print("V1:\n", V1)
    print("V2:\n", V2)

    # Inner product: Use np.inner, .dot, .matmul for 1D list, ndarray
    print("")
    print("Inner product")
    inner = np.inner(V1, V2)
    print("np.inner V1 dot V2 = {}".format(inner))
    dot = np.dot(V1, V2)
    print("np.dot V1 dot V2 = {}".format(dot))

    # Outer product
    V3 = np.cross(V1, V2)
    print("")
    print("Outer product (DO NOT use np.outer()")
    print("np.cross V1 x V2:")
    print(V3)

    # Inverse matrix: linalg.inv for list
    A_i = npalg.inv(A)
    check = A @ A_i
    print("")
    print("np.inv A^-1:")
    print(A_i)
    print("A*A.i:")
    print(check)

    # Determinant
    arr = np.array([[0, 1], [2, 3]])
    det = np.linalg.det(arr)
    print("")
    print("determinant for A = {}".format(det))

    # Eigen problem
    lA, vA = npalg.eig(A)
    print("")
    print("nympy.linalg.eig(A):")
    print("Eigen values: ", lA)
    print("Eigen vectors:")
    print(vA)

    # Solve simultaneous linear equations
    print("")
    print("Solve simultaneous linear equations AX = B:")
    X = np.linalg.solve(A, B)
    print(X)

    # LU decomposition
    print("")
    print("LU decomposition A = PLU:")
    print("A:\n", A)
    P, L, U = spalg.lu(A)
    check_lu = P @ L @ U
    print("P:\n", P)
    print("L:\n", L)
    print("U:\n", U)
    print("P@L@U:\n", check_lu)

    # Cholesky decomosition: S = LL^T
    print("Cholesky decompositioin S = LL^T:")
    L_chol = npalg.cholesky(S)
    check_chol = L_chol @ L_chol.T
    print("")
    print("S:\n", S)
    print("L:\n", L_chol)
    print("L*L^T:\n", check_chol)

    # QR decomposition: A = QR
    # good for anomarly matrix
    Q, R = spalg.qr(A)
    check_qr = Q @ R
    print("")
    print("QR decomposition:")
    print("A:\n", A)
    print("Q:\n", Q)
    print("R:\n", R)
    print("Q@R:\n", check_qr)

    print("")


if __name__ == "__main__":
    main()