工学のための線形代数学

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt

# ガウス消去法による線形方程式の解法
def gaussian_elimination(A, b):
    """
    ガウス消去法を用いて線形方程式 Ax = b を解く。
    A: 係数行列
    b: 右辺ベクトル
    """
    A = np.array(A, dtype=float)
    b = np.array(b, dtype=float)
    n = len(b)
    # 前進消去
    for i in range(n):
        # 対角要素を1にする
        A[i] = A[i] / A[i, i]
        b[i] = b[i] / A[i, i]
        for j in range(i + 1, n):
            factor = A[j, i]
            A[j, i:] -= factor * A[i, i:]
            b[j] -= factor * b[i]
    # 後退代入
    x = np.zeros(n)
    for i in reversed(range(n)):
        x[i] = b[i] - np.dot(A[i, i + 1:], x[i + 1:])
    return x

# LU分解
def lu_decomposition(A):
    """
    LU分解を行う。
    A: 係数行列
    """
    P, L, U = la.lu(A)
    return L, U, P

def solve_with_lu(L, U, P, b):
    """
    LU分解を用いて線形方程式 Ax = b を解く。
    L: 下三角行列
    U: 上三角行列
    P: ピボット行列
    b: 右辺ベクトル
    """
    Pb = np.dot(P, b)
    y = np.linalg.solve(L, Pb)
    x = np.linalg.solve(U, y)
    return x

# 例: 線形方程式の系を解く
A = np.array([[2, -1, 1], [1, 3, 2], [1, 1, 2]])
b = np.array([2, 5, 3])

# ガウス消去法による解
x_gaussian = gaussian_elimination(A.copy(), b.copy())

# LU分解
L, U, P = lu_decomposition(A)

# LU分解による解
x_lu = solve_with_lu(L, U, P, b)

# 結果の出力
print("ガウス消去法による解:", x_gaussian)
print("LU分解による解:", x_lu)

# 結果を比較するためのプロット
fig, ax = plt.subplots(1, 2, figsize=(14, 6))

# ガウス消去法の結果をプロット
ax[0].bar(range(len(x_gaussian)), x_gaussian, color='blue', alpha=0.7)
ax[0].set_title('Gaussian Elimination Results')
ax[0].set_xlabel('Variable Index')
ax[0].set_ylabel('Value')

# LU分解の結果をプロット
ax[1].bar(range(len(x_lu)), x_lu, color='green', alpha=0.7)
ax[1].set_title('LU Decomposition Results')
ax[1].set_xlabel('Variable Index')
ax[1].set_ylabel('Value')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# ベクトル空間と部分空間の定義
def basis_and_dimension(vectors):
    """
    与えられたベクトルのリストから基底と次元を計算する。
    vectors: ベクトルのリスト (numpy配列)
    """
    matrix = np.array(vectors).T
    rank = np.linalg.matrix_rank(matrix)
    U, S, Vt = np.linalg.svd(matrix, full_matrices=False)
    basis = Vt[:rank]
    return basis.T, rank

# 線形変換の定義と行列による表現
def linear_transformation(matrix, vectors):
    """
    線形変換を行う。
    matrix: 変換行列
    vectors: 入力ベクトルのリスト (numpy配列)
    """
    return [np.dot(matrix, v) for v in vectors]

# 固有値と固有ベクトルの計算
def eigen_decomposition(matrix):
    """
    行列の固有値と固有ベクトルを計算する。
    matrix: 対角化する行列
    """
    eigenvalues, eigenvectors = np.linalg.eig(matrix)
    return eigenvalues, eigenvectors

# 行列の対角化
def diagonalize_matrix(matrix):
    """
    行列を対角化する。
    matrix: 対角化する行列
    """
    eigenvalues, eigenvectors = eigen_decomposition(matrix)
    D = np.diag(eigenvalues)
    P = eigenvectors
    P_inv = np.linalg.inv(P)
    A_diagonalized = P @ D @ P_inv
    return A_diagonalized

# 例: ベクトル空間
vectors = [[1, 0], [0, 1], [1, 1]]
basis, dimension = basis_and_dimension(vectors)

# 例: 線形変換
transformation_matrix = np.array([[2, 0], [0, 3]])
transformed_vectors = linear_transformation(transformation_matrix, vectors)

# 例: 固有値と固有ベクトル
matrix = np.array([[4, 1], [2, 3]])
eigenvalues, eigenvectors = eigen_decomposition(matrix)
A_diagonalized = diagonalize_matrix(matrix)

# 結果の出力
print("基底:\n", basis)
print("次元:", dimension)
print("変換後のベクトル:\n", transformed_vectors)
print("固有値:", eigenvalues)
print("固有ベクトル:\n", eigenvectors)
print("対角化された行列:\n", A_diagonalized)

# 結果のプロット
fig, ax = plt.subplots(1, 2, figsize=(14, 6))

# 基底ベクトルのプロット
ax[0].quiver(0, 0, basis[0, 0], basis[1, 0], angles='xy', scale_units='xy', scale=1, color='r', label='Basis Vector 1')
ax[0].quiver(0, 0, basis[0, 1], basis[1, 1], angles='xy', scale_units='xy', scale=1, color='b', label='Basis Vector 2')
ax[0].set_xlim(-2, 2)
ax[0].set_ylim(-2, 2)
ax[0].set_title('Basis Vectors')
ax[0].set_xlabel('x')
ax[0].set_ylabel('y')
ax[0].legend()
ax[0].grid()

# 固有ベクトルのプロット
ax[1].quiver(0, 0, eigenvectors[0, 0], eigenvectors[1, 0], angles='xy', scale_units='xy', scale=1, color='g', label='Eigenvector 1')
ax[1].quiver(0, 0, eigenvectors[0, 1], eigenvectors[1, 1], angles='xy', scale_units='xy', scale=1, color='m', label='Eigenvector 2')
ax[1].set_xlim(-2, 2)
ax[1].set_ylim(-2, 2)
ax[1].set_title('Eigenvectors')
ax[1].set_xlabel('x')
ax[1].set_ylabel('y')
ax[1].legend()
ax[1].grid()

plt.tight_layout()
plt.show()

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt

# 固有値と固有ベクトルの計算
def eigen_decomposition(matrix):
    """
    行列の固有値と固有ベクトルを計算する。
    matrix: 対角化する行列
    """
    eigenvalues, eigenvectors = np.linalg.eig(matrix)
    return eigenvalues, eigenvectors

# 行列の対角化
def diagonalize_matrix(matrix):
    """
    行列を対角化する。
    matrix: 対角化する行列
    """
    eigenvalues, eigenvectors = eigen_decomposition(matrix)
    D = np.diag(eigenvalues)
    P = eigenvectors
    P_inv = np.linalg.inv(P)
    A_diagonalized = P @ D @ P_inv
    return A_diagonalized

# ジョルダン標準形の計算
def jordan_normal_form(matrix):
    """
    行列のジョルダン標準形を計算する。
    matrix: ジョルダン標準形を求める行列
    """
    J, P = la.jordan(matrix)
    return J, P

# 例: 固有値と固有ベクトルの計算
matrix = np.array([[1, 2], [3, 4]])
eigenvalues, eigenvectors = eigen_decomposition(matrix)
A_diagonalized = diagonalize_matrix(matrix)

# 例: ジョルダン標準形の計算
jordan_matrix, jordan_transform = jordan_normal_form(matrix)

# 結果の出力
print("固有値:", eigenvalues)
print("固有ベクトル:\n", eigenvectors)
print("対角化された行列:\n", A_diagonalized)
print("ジョルダン標準形:\n", jordan_matrix)
print("ジョルダン変換行列:\n", jordan_transform)

# 結果のプロット（固有ベクトルのプロット）
fig, ax = plt.subplots(figsize=(7, 7))

# 固有ベクトルのプロット
ax.quiver(0, 0, eigenvectors[0, 0], eigenvectors[1, 0], angles='xy', scale_units='xy', scale=1, color='g', label='Eigenvector 1')
ax.quiver(0, 0, eigenvectors[0, 1], eigenvectors[1, 1], angles='xy', scale_units='xy', scale=1, color='b', label='Eigenvector 2')
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_title('Eigenvectors')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.legend()
ax.grid()

plt.show()

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt

# 4.1 LU分解
def lu_decomposition(matrix):
    """
    LU分解を計算する。
    matrix: 分解する行列
    """
    P, L, U = la.lu(matrix)
    return P, L, U

# 4.2 QR分解
def qr_decomposition(matrix):
    """
    QR分解を計算する。
    matrix: 分解する行列
    """
    Q, R = np.linalg.qr(matrix)
    return Q, R

def least_squares_solution(A, b):
    """
    最小二乗法による解を計算する。
    A: 行列
    b: ベクトル
    """
    Q, R = qr_decomposition(A)
    y = Q.T @ b
    x = np.linalg.solve(R, y)
    return x

# 4.3 特異値分解 (SVD)
def svd_decomposition(matrix):
    """
    特異値分解を計算する。
    matrix: 分解する行列
    """
    U, S, Vt = np.linalg.svd(matrix, full_matrices=False)
    return U, S, Vt

def data_compression(matrix, k):
    """
    特異値分解を用いたデータ圧縮を行う。
    matrix: 圧縮する行列
    k: 次元数
    """
    U, S, Vt = svd_decomposition(matrix)
    S_k = np.diag(S[:k])
    U_k = U[:, :k]
    Vt_k = Vt[:k, :]
    compressed_matrix = U_k @ S_k @ Vt_k
    return compressed_matrix

# 例: LU分解
matrix = np.array([[4, 2], [3, 1]])
P, L, U = lu_decomposition(matrix)

# 例: QR分解と最小二乗法
A = np.array([[1, 2], [3, 4], [5, 6]])
b = np.array([7, 8, 9])
x_ls = least_squares_solution(A, b)

# 例: 特異値分解とデータ圧縮
data_matrix = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
k = 2
compressed_data_matrix = data_compression(data_matrix, k)

# 結果の出力
print("LU分解:")
print("P:\n", P)
print("L:\n", L)
print("U:\n", U)

print("\nQR分解:")
print("Q:\n", qr_decomposition(A)[0])
print("R:\n", qr_decomposition(A)[1])

print("\n最小二乗法の解:")
print("x:\n", x_ls)

print("\n特異値分解:")
U_svd, S_svd, Vt_svd = svd_decomposition(data_matrix)
print("U:\n", U_svd)
print("S:\n", S_svd)
print("Vt:\n", Vt_svd)

print("\nデータ圧縮:")
print("圧縮された行列:\n", compressed_data_matrix)

# 結果のプロット（SVDの特異値）
plt.figure(figsize=(8, 6))
plt.plot(S_svd, 'o-', label='Singular Values')
plt.xlabel('Index')
plt.ylabel('Singular Value')
plt.title('Singular Values from SVD')
plt.legend()
plt.grid()
plt.show()

import numpy as np
from scipy.optimize import linprog
import matplotlib.pyplot as plt

# 5.1 最適化の基礎

# 線形最適化問題の定式化
# 目的関数: c^T * x の最小化
# 制約条件: A * x <= b

# 目的関数の係数
c = np.array([-1, -2])  # 例: 最大化問題のため、係数は負にする

# 制約条件の係数行列
A = np.array([[2, 1], [1, 2]])

# 制約条件の右辺
b = np.array([20, 30])

# 制約条件の等式/不等式の種類（全て不等式）
# '<=', '>=', '=' で指定
# ここでは全て '<=' です
constraints = 'ineq'

# 線形計画問題の解法
result = linprog(c, A_ub=A, b_ub=b, method='simplex')

# 結果の表示
print("最適化結果:")
print("最適解:", result.x)
print("最適値:", result.fun)
print("ステータス:", result.message)

# 制約条件と目的関数のプロット
# 可視化用のデータ生成
x = np.linspace(0, 10, 400)
y1 = (20 - 2 * x) / 1
y2 = (30 - 1 * x) / 2

plt.figure(figsize=(8, 6))
plt.plot(x, y1, label=r'$2x_1 + x_2 \leq 20$')
plt.plot(x, y2, label=r'$x_1 + 2x_2 \leq 30$')
plt.fill_between(x, np.minimum(y1, y2), 0, alpha=0.2, color='gray', label='Feasible Region')

plt.scatter(result.x[0], result.x[1], color='red', label='Optimal Solution')
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.title('Linear Programming Feasible Region and Optimal Solution')
plt.legend()
plt.grid()
plt.show()

import numpy as np
from scipy.optimize import linprog
import matplotlib.pyplot as plt

# 5.2 線形計画法
# シンプレックス法
# 双対性理論

# プライマリ問題の定義
c = np.array([-1, -2])  # 目的関数係数（最大化のため負に設定）
A = np.array([[2, 1], [1, 2]])
b = np.array([20, 30])

# 双対問題の定義
c_dual = np.array([20, 30])  # 双対の目的関数係数
A_dual = np.array([[2, 1], [1, 2]])
b_dual = np.array([-1, -2])  # 双対の制約条件の右辺（プライマリの係数）

# プライマリ問題を解く
result_primal = linprog(c, A_ub=A, b_ub=b, method='simplex')
print("プライマリ問題の最適解:")
print("最適解:", result_primal.x)
print("最適値:", result_primal.fun)

# 双対問題を解く
result_dual = linprog(c_dual, A_ub=-A_dual, b_ub=-b_dual, method='simplex')
print("双対問題の最適解:")
print("最適解:", result_dual.x)
print("最適値:", result_dual.fun)

# プライマリ問題と双対問題の可視化
x = np.linspace(0, 10, 400)
y1 = (20 - 2 * x) / 1
y2 = (30 - 1 * x) / 2

plt.figure(figsize=(12, 8))
plt.subplot(1, 2, 1)
plt.plot(x, y1, label=r'$2x_1 + x_2 \leq 20$')
plt.plot(x, y2, label=r'$x_1 + 2x_2 \leq 30$')
plt.fill_between(x, np.minimum(y1, y2), 0, alpha=0.2, color='gray', label='Feasible Region')
plt.scatter(result_primal.x[0], result_primal.x[1], color='red', label='Primal Optimal Solution')
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.title('Primal Problem')
plt.legend()
plt.grid()

plt.subplot(1, 2, 2)
plt.plot(x, -(-20 + 2 * x) / 1, label=r'$2x_1 + x_2 \geq -20$')
plt.plot(x, -(-30 + 1 * x) / 2, label=r'$x_1 + 2x_2 \geq -30$')
plt.fill_between(x, np.maximum(-(-20 + 2 * x) / 1, -(-30 + 1 * x) / 2), 0, alpha=0.2, color='lightblue', label='Feasible Region')
plt.scatter(result_dual.x[0], result_dual.x[1], color='green', label='Dual Optimal Solution')
plt.xlabel(r'$y_1$')
plt.ylabel(r'$y_2$')
plt.title('Dual Problem')
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt

# 数値解法の精度と行列計算の効率化

# 行列Aとベクトルb
A = np.array([[4, 1], [3, 2]])
b = np.array([1, 2])

# ガウス消去法で解を求める
x_gauss = np.linalg.solve(A, b)
print("Gaussian Elimination Solution:", x_gauss)

# LU分解
P, L, U = la.lu(A)
x_lu = la.solve(U, la.solve(L, np.dot(P.T, b)))
print("LU Decomposition Solution:", x_lu)

# QR分解
Q, R = np.linalg.qr(A)
x_qr = np.linalg.solve(R, np.dot(Q.T, b))
print("QR Decomposition Solution:", x_qr)

# 精度比較
def compute_residual(A, x, b):
    return np.linalg.norm(np.dot(A, x) - b)

residual_gauss = compute_residual(A, x_gauss, b)
residual_lu = compute_residual(A, x_lu, b)
residual_qr = compute_residual(A, x_qr, b)

print("Residual (Gaussian Elimination):", residual_gauss)
print("Residual (LU Decomposition):", residual_lu)
print("Residual (QR Decomposition):", residual_qr)

# 行列計算の効率化のプロット
matrix_sizes = np.arange(10, 510, 50)
lu_times = []
qr_times = []

for size in matrix_sizes:
    A_large = np.random.rand(size, size)
    b_large = np.random.rand(size)

    # LU分解の計算時間
    start_time = np.time()
    la.lu(A_large)
    lu_times.append(np.time() - start_time)

    # QR分解の計算時間
    start_time = np.time()
    np.linalg.qr(A_large)
    qr_times.append(np.time() - start_time)

plt.figure(figsize=(12, 6))
plt.plot(matrix_sizes, lu_times, label='LU Decomposition Time', marker='o')
plt.plot(matrix_sizes, qr_times, label='QR Decomposition Time', marker='o')
plt.xlabel('Matrix Size')
plt.ylabel('Time (seconds)')
plt.title('Performance of Matrix Decompositions')
plt.legend()
plt.grid()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 線形変換の定義
def linear_transformation(A, x):
    return np.dot(A, x)

# 行列の定義
A = np.array([[2, 0], [0, 3]])
B = np.array([[1, 2], [3, 4]])

# ベクトルの定義
x = np.array([1, 2])

# 線形変換
T_x = linear_transformation(A, x)
S_x = linear_transformation(B, x)

# 合成変換
T_S_x = linear_transformation(A, S_x)

# 逆変換
A_inv = np.linalg.inv(A)
B_inv = np.linalg.inv(B)
T_inv_x = linear_transformation(A_inv, x)
S_inv_x = linear_transformation(B_inv, x)

# プロット設定
plt.figure(figsize=(12, 8))

# ベクトルのプロット
plt.subplot(2, 2, 1)
plt.quiver(0, 0, x[0], x[1], angles='xy', scale_units='xy', scale=1, color='r')
plt.title('Original Vector')
plt.xlim(-1, 10)
plt.ylim(-1, 10)
plt.grid()

# 線形変換後のベクトル
plt.subplot(2, 2, 2)
plt.quiver(0, 0, T_x[0], T_x[1], angles='xy', scale_units='xy', scale=1, color='b', label='A(x)')
plt.quiver(0, 0, S_x[0], S_x[1], angles='xy', scale_units='xy', scale=1, color='g', label='B(x)')
plt.title('Transformed Vectors')
plt.xlim(-1, 10)
plt.ylim(-1, 10)
plt.legend()
plt.grid()

# 合成変換後のベクトル
plt.subplot(2, 2, 3)
plt.quiver(0, 0, T_S_x[0], T_S_x[1], angles='xy', scale_units='xy', scale=1, color='m', label='T(S(x))')
plt.title('Composite Transformation')
plt.xlim(-1, 10)
plt.ylim(-1, 10)
plt.legend()
plt.grid()

# 逆変換後のベクトル
plt.subplot(2, 2, 4)
plt.quiver(0, 0, T_inv_x[0], T_inv_x[1], angles='xy', scale_units='xy', scale=1, color='c', label='A^(-1)(x)')
plt.quiver(0, 0, S_inv_x[0], S_inv_x[1], angles='xy', scale_units='xy', scale=1, color='y', label='B^(-1)(x)')
plt.title('Inverse Transformations')
plt.xlim(-1, 10)
plt.ylim(-1, 10)
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

import numpy as np
import scipy.linalg as linalg
import matplotlib.pyplot as plt

# サンプル行列の定義
A = np.array([[4, 2], [3, 6]])
B = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

# LU分解
P, L, U = linalg.lu(A)

# QR分解
Q, R = linalg.qr(B)

# 特異値分解 (SVD)
U_svd, S, Vt = linalg.svd(B)

# プロット設定
fig, axs = plt.subplots(2, 2, figsize=(12, 10))

# LU分解の表示
axs[0, 0].imshow(L, cmap='viridis', interpolation='none')
axs[0, 0].set_title('L Matrix (LU Decomposition)')
axs[0, 0].grid(False)

axs[0, 1].imshow(U, cmap='viridis', interpolation='none')
axs[0, 1].set_title('U Matrix (LU Decomposition)')
axs[0, 1].grid(False)

# QR分解の表示
axs[1, 0].imshow(Q, cmap='viridis', interpolation='none')
axs[1, 0].set_title('Q Matrix (QR Decomposition)')
axs[1, 0].grid(False)

axs[1, 1].imshow(R, cmap='viridis', interpolation='none')
axs[1, 1].set_title('R Matrix (QR Decomposition)')
axs[1, 1].grid(False)

plt.tight_layout()
plt.show()

# 特異値分解の表示
fig, axs = plt.subplots(1, 2, figsize=(12, 5))

axs[0].imshow(U_svd, cmap='viridis', interpolation='none')
axs[0].set_title('U Matrix (SVD)')
axs[0].grid(False)

axs[1].imshow(Vt, cmap='viridis', interpolation='none')
axs[1].set_title('V^T Matrix (SVD)')
axs[1].grid(False)

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# サンプルデータの生成
np.random.seed(0)
X = np.random.rand(100, 5)  # 100サンプル、5次元のデータ

# データの標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# PCAの実行
pca = PCA(n_components=2)  # 2次元に削減
X_pca = pca.fit_transform(X_scaled)

# 元のデータとPCA後のデータのプロット
fig, axs = plt.subplots(1, 2, figsize=(12, 5))

# 元のデータの散布図
axs[0].scatter(X_scaled[:, 0], X_scaled[:, 1], alpha=0.7)
axs[0].set_title('Original Data (Scaled)')
axs[0].set_xlabel('Feature 1')
axs[0].set_ylabel('Feature 2')

# PCA後のデータの散布図
axs[1].scatter(X_pca[:, 0], X_pca[:, 1], alpha=0.7)
axs[1].set_title('PCA Reduced Data')
axs[1].set_xlabel('Principal Component 1')
axs[1].set_ylabel('Principal Component 2')

plt.tight_layout()
plt.show()

# PCAの主成分と寄与率の表示
print("Principal Components:")
print(pca.components_)

print("\nExplained Variance Ratio:")
print(pca.explained_variance_ratio_)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import classification_report
from sklearn.svm import SVC

# データセットの読み込み
data = load_iris()
X = data.data
y = data.target
feature_names = data.feature_names

# データの標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# PCAの実行
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

# LDAの実行
lda = LDA(n_components=2)
X_lda = lda.fit_transform(X_scaled, y)

# データの分割
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, test_size=0.3, random_state=42)

# PCA後のデータでモデル訓練と評価
clf_pca = SVC(kernel='linear')
clf_pca.fit(X_train, y_train)
y_pred_pca = clf_pca.predict(X_test)

# LDA後のデータでモデル訓練と評価
clf_lda = SVC(kernel='linear')
clf_lda.fit(X_train, y_train)
y_pred_lda = clf_lda.predict(X_test)

# PCAとLDAのプロット
fig, axs = plt.subplots(1, 2, figsize=(14, 6))

# PCAのプロット
scatter_pca = axs[0].scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap='viridis')
axs[0].set_title('PCA Reduced Data')
axs[0].set_xlabel('Principal Component 1')
axs[0].set_ylabel('Principal Component 2')
legend1 = axs[0].legend(*scatter_pca.legend_elements(), title="Classes")
axs[0].add_artist(legend1)

# LDAのプロット
scatter_lda = axs[1].scatter(X_lda[:, 0], X_lda[:, 1], c=y, cmap='viridis')
axs[1].set_title('LDA Reduced Data')
axs[1].set_xlabel('LD 1')
axs[1].set_ylabel('LD 2')
legend2 = axs[1].legend(*scatter_lda.legend_elements(), title="Classes")
axs[1].add_artist(legend2)

plt.tight_layout()
plt.show()

# モデル評価
print("PCA + SVM Classification Report:")
print(classification_report(y_test, y_pred_pca))

print("\nLDA + SVM Classification Report:")
print(classification_report(y_test, y_pred_lda))

import numpy as np
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# サンプルデータの生成
np.random.seed(0)
data = np.random.rand(100, 5)  # 100サンプル、5次元のデータ

# PCAの実行
pca = PCA(n_components=2)  # 2次元に削減
transformed_data = pca.fit_transform(data)

# 固有値と固有ベクトルの取得
explained_variance = pca.explained_variance_ratio_
components = pca.components_

# プロット
plt.figure(figsize=(8, 6))
plt.scatter(transformed_data[:, 0], transformed_data[:, 1], alpha=0.7)
plt.title('PCA Result')
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.grid(True)
plt.show()

# 固有値と固有ベクトルの表示
print("Explained Variance Ratio (Eigenvalues):", explained_variance)
print("Principal Components (Eigenvectors):\n", components)

import numpy as np
import matplotlib.pyplot as plt

# Example adjacency matrix of a gene network (undirected graph)
adj_matrix = np.array([
    [0, 1, 0, 1],
    [1, 0, 1, 0],
    [0, 1, 0, 1],
    [1, 0, 1, 0]
])

# Calculate eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(adj_matrix)

# Plot eigenvalues
plt.figure(figsize=(8, 6))
plt.scatter(range(len(eigenvalues)), eigenvalues, color='blue', marker='o')
plt.title('Eigenvalues of the Gene Network Adjacency Matrix')
plt.xlabel('Index')
plt.ylabel('Eigenvalue')
plt.grid(True)
plt.show()

# Display eigenvalues and eigenvectors
print("Eigenvalues:\n", eigenvalues)
print("Eigenvectors:\n", eigenvectors)

import numpy as np
import matplotlib.pyplot as plt

# Example Jacobian matrix of an enzyme reaction kinetics model
jacobian_matrix = np.array([
    [-0.1, 0.05],
    [0.02, -0.07]
])

# Calculate eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(jacobian_matrix)

# Plot eigenvalues
plt.figure(figsize=(8, 6))
plt.scatter(range(len(eigenvalues)), eigenvalues, color='red', marker='o')
plt.title('Eigenvalues of the Jacobian Matrix for Enzyme Reaction')
plt.xlabel('Index')
plt.ylabel('Eigenvalue')
plt.grid(True)
plt.show()

# Display eigenvalues and eigenvectors
print("Eigenvalues:\n", eigenvalues)
print("Eigenvectors:\n", eigenvectors)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# データの生成（例として2次元データを生成）
np.random.seed(0)
X = np.random.randn(100, 2)

# データの標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# PCAの実施
pca = PCA(n_components=1)
X_pca = pca.fit_transform(X_scaled)

# 結果のプロット
plt.figure(figsize=(8, 6))
plt.scatter(X_scaled[:, 0], X_scaled[:, 1], label='Original Data')
plt.scatter(X_pca, np.zeros_like(X_pca), label='PCA Reduced Data', color='red')
plt.xlabel('Feature 1')
plt.ylabel('Feature 2')
plt.title('PCA Dimension Reduction')
plt.legend()
plt.grid(True)
plt.show()

print('Explained variance ratio:', pca.explained_variance_ratio_)
print('Principal Components:\n', pca.components_)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 線形システムの定義
def system(x, t):
    A = np.array([[-1, 2], [-2, -1]])
    return A @ x

# リアプノフ関数の定義
def lyapunov_function(x):
    return x.T @ P @ x

# リアプノフ関数の時間変化
def lyapunov_derivative(x):
    return x.T @ (A.T @ P + P @ A) @ x

# システムとリアプノフ関数の設定
A = np.array([[-1, 2], [-2, -1]])
P = np.linalg.solve(A.T @ A, np.eye(2))  # P を適切に選択
x0 = np.array([1, 0])  # 初期状態

# 時間の設定
t = np.linspace(0, 10, 100)

# システムの解の計算
x = odeint(system, x0, t)

# リアプノフ関数の値の計算
lyapunov_values = np.array([lyapunov_function(xi) for xi in x])

# プロット
plt.figure(figsize=(10, 6))
plt.plot(t, lyapunov_values, label='Lyapunov Function V(x)')
plt.xlabel('Time')
plt.ylabel('Lyapunov Function Value')
plt.title('Lyapunov Function Over Time')
plt.legend()
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# シミュレーションパラメータ
num_particles = 100
num_steps = 1000
dt = 0.01

# 初期位置と速度
positions = np.random.rand(num_particles, 3)
velocities = np.random.randn(num_particles, 3)

# 力の計算（シンプルな例）
def compute_forces(positions):
    forces = np.zeros_like(positions)
    for i in range(num_particles):
        for j in range(i + 1, num_particles):
            r = positions[j] - positions[i]
            r_norm = np.linalg.norm(r)
            if r_norm > 0:
                force = r / r_norm ** 3
                forces[i] += force
                forces[j] -= force
    return forces

# シミュレーションループ
trajectory = []
for step in range(num_steps):
    forces = compute_forces(positions)
    velocities += forces * dt
    positions += velocities * dt
    trajectory.append(positions.copy())

# 結果のプロット
trajectory = np.array(trajectory)
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')

for i in range(num_particles):
    ax.plot(trajectory[:, i, 0], trajectory[:, i, 1], trajectory[:, i, 2])

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.title('Molecular Dynamics Simulation')
plt.show()

import numpy as np

def partial_pivoting_gaussian_elimination(A, b):
    n = len(b)
    
    # 拡大係数行列の作成
    Ab = np.hstack([A, b.reshape(-1, 1)])
    
    for i in range(n):
        # ピボット選択
        max_row = np.argmax(np.abs(Ab[i:n, i])) + i
        if i != max_row:
            # 行の交換
            Ab[[i, max_row]] = Ab[[max_row, i]]
        
        # ガウス消去法
        for j in range(i + 1, n):
            ratio = Ab[j, i] / Ab[i, i]
            Ab[j, i:] -= ratio * Ab[i, i:]
    
    # 後退代入法
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - np.dot(Ab[i, i + 1:n], x[i + 1:n])) / Ab[i, i]
    
    return x

# 例
A = np.array([[2, -1, 1], [3, 3, 9], [3, 3, 5]], dtype=float)
b = np.array([1, 0, 4], dtype=float)

x = partial_pivoting_gaussian_elimination(A, b)
print("Solution x:", x)

import numpy as np
import matplotlib.pyplot as plt

# データの生成
np.random.seed(0)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

# 行列 A とベクトル b の作成
A = np.hstack([np.ones((100, 1)), X])  # バイアス項の追加
b = y

# 正規方程式の計算
AtA = A.T @ A
Atb = A.T @ b
x = np.linalg.solve(AtA, Atb)

# 予測値の計算
y_pred = A @ x

# プロット
plt.scatter(X, y, label='Data points')
plt.plot(X, y_pred, color='red', label='Regression line')
plt.xlabel('X')
plt.ylabel('y')
plt.legend()
plt.title('Linear Regression using Normal Equations')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error

# データの生成
np.random.seed(0)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X + np.random.randn(100, 1)

# バイアス項の追加
X_b = np.c_[np.ones((100, 1)), X]

# 正規方程式の計算
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)

# 予測値の計算
X_new = np.array([[0], [2]])
X_new_b = np.c_[np.ones((2, 1)), X_new]
y_predict = X_new_b.dot(theta_best)

# プロット
plt.plot(X_new, y_predict, "r-", label="Prediction")
plt.plot(X, y, "b.", label="Data points")
plt.xlabel("X")
plt.ylabel("y")
plt.legend()
plt.title("Linear Regression using Normal Equations")
plt.show()

# モデルの評価
y_pred = X_b.dot(theta_best)
mse = mean_squared_error(y, y_pred)
print("Mean Squared Error:", mse)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import precision_score, recall_score, f1_score, accuracy_score
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression

# データの生成
X, y = make_classification(n_samples=1000, n_features=20, n_classes=2, random_state=42)

# トレーニングとテストデータに分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# ロジスティック回帰モデルのトレーニング
model = LogisticRegression()
model.fit(X_train, y_train)

# 予測
y_pred = model.predict(X_test)

# 評価指標の計算
accuracy = accuracy_score(y_test, y_pred)
precision = precision_score(y_test, y_pred)
recall = recall_score(y_test, y_pred)
f1 = f1_score(y_test, y_pred)

# 評価指標を辞書に格納
metrics = {
    'Accuracy': accuracy,
    'Precision': precision,
    'Recall': recall,
    'F1 Score': f1
}

# プロット
plt.figure(figsize=(8, 6))
plt.bar(metrics.keys(), metrics.values(), color=['blue', 'green', 'red', 'purple'])
plt.ylim(0, 1)
plt.title('Model Evaluation Metrics')
plt.ylabel('Score')
plt.xlabel('Metrics')
plt.grid(axis='y')

# 各バーにスコアを表示
for i, (metric, score) in enumerate(metrics.items()):
    plt.text(i, score + 0.02, f'{score:.2f}', ha='center')

plt.show()

# 必要なライブラリのインポート
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

# データの生成
np.random.seed(0)
mean = [0, 0]
cov = [[3, 1], [1, 2]]  # 分散共分散行列
X = np.random.multivariate_normal(mean, cov, 100)

# データのプロット
plt.scatter(X[:, 0], X[:, 1], alpha=0.5)
plt.title('Original Data')
plt.xlabel('X1')
plt.ylabel('X2')
plt.axis('equal')
plt.show()

# 分散共分散行列の計算
cov_matrix = np.cov(X.T)
print("Covariance Matrix:\n", cov_matrix)

# 固有値と固有ベクトルの計算
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
print("Eigenvalues:\n", eigenvalues)
print("Eigenvectors:\n", eigenvectors)

# PCAの適用
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)

# 主成分軸をプロット
plt.scatter(X[:, 0], X[:, 1], alpha=0.5, label='Original Data')
for length, vector in zip(eigenvalues, eigenvectors.T):
    v = vector * 3 * np.sqrt(length)
    plt.quiver(*mean, *v, angles='xy', scale_units='xy', scale=1, color='r', alpha=0.8)

plt.title('Principal Components')
plt.xlabel('X1')
plt.ylabel('X2')
plt.axis('equal')
plt.legend()
plt.show()

# PCA後のデータプロット
plt.scatter(X_pca[:, 0], X_pca[:, 1], alpha=0.5)
plt.title('Data after PCA')
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.axis('equal')
plt.show()

# 必要なライブラリのインポート
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import TruncatedSVD
from sklearn.preprocessing import StandardScaler

# データの生成
np.random.seed(0)
data = np.random.randn(100, 50)  # 100サンプル、50特徴

# 標準化
scaler = StandardScaler()
data_standardized = scaler.fit_transform(data)

# 特異値分解の適用
svd = TruncatedSVD(n_components=10)
data_svd = svd.fit_transform(data_standardized)

# 特異値の取得
singular_values = svd.singular_values_
print("Singular Values:\n", singular_values)

# 元のデータの可視化
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.title("Original Data")
plt.imshow(data, aspect='auto', cmap='viridis')
plt.colorbar()
plt.xlabel("Features")
plt.ylabel("Samples")

# SVD後のデータの可視化
plt.subplot(1, 2, 2)
plt.title("Data after SVD")
plt.imshow(data_svd, aspect='auto', cmap='viridis')
plt.colorbar()
plt.xlabel("Reduced Features")
plt.ylabel("Samples")
plt.tight_layout()
plt.show()

# 深層学習における特異値の役割の例
def initialize_weights(shape):
    return np.random.randn(*shape) * np.sqrt(2 / sum(shape))

# 勾配消失・爆発の例
layer_sizes = [50, 100, 200, 100, 50]
weights = [initialize_weights((layer_sizes[i], layer_sizes[i+1])) for i in range(len(layer_sizes) - 1)]

# 重みの積の特異値の計算
product_matrix = weights[0]
for W in weights[1:]:
    product_matrix = np.dot(product_matrix, W)

_, singular_values, _ = np.linalg.svd(product_matrix)
print("Singular Values of the Weight Product Matrix:\n", singular_values)

# 特異値のプロット
plt.figure(figsize=(8, 6))
plt.plot(singular_values, 'o-')
plt.title("Singular Values of the Weight Product Matrix")
plt.xlabel("Index")
plt.ylabel("Singular Value")
plt.yscale('log')
plt.show()

# 必要なライブラリのインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
import torch
import torch.nn as nn
import torch.optim as optim

# Chapter 18: Regression Analysis and Pseudoinverse
print("Chapter 18: Regression Analysis and Pseudoinverse")

# データの生成
np.random.seed(0)
X = 2 * np.random.rand(100, 2)
y = 4 + 3 * X[:, 0] + 5 * X[:, 1] + np.random.randn(100)

# 拡張行列
X_b = np.c_[np.ones((100, 1)), X]

# 擬似逆行列を用いた回帰分析
theta_best = np.linalg.pinv(X_b).dot(y)
print("Estimated coefficients:", theta_best)

# 結果のプロット
plt.figure(figsize=(8, 6))
plt.scatter(X[:, 0], y, label='Data points')
plt.plot(X[:, 0], X_b.dot(theta_best), color='red', label='Regression line')
plt.xlabel("X1")
plt.ylabel("y")
plt.title("Regression Analysis with Pseudoinverse")
plt.legend()
plt.show()

# Chapter 19: Multivariate Normal Distribution and Its Integration
print("Chapter 19: Multivariate Normal Distribution and Its Integration")

# 多変量正規分布の定義
mean = np.array([0, 0])
cov = np.array([[1, 0.8], [0.8, 1]])
rv = multivariate_normal(mean, cov)

# 分布のサンプル生成とプロット
X = np.random.multivariate_normal(mean, cov, 500)
plt.figure(figsize=(8, 6))
plt.scatter(X[:, 0], X[:, 1], alpha=0.5, label='Samples')
plt.title("Multivariate Normal Distribution")
plt.xlabel("X1")
plt.ylabel("X2")
plt.axis('equal')
plt.legend()
plt.show()

# Chapter 20: Generative Models and Variational Free Energy
print("Chapter 20: Generative Models and Variational Free Energy")

# 簡単な生成モデルの実装 (Variational Autoencoder)
class VAE(nn.Module):
    def __init__(self, input_dim, hidden_dim, latent_dim):
        super(VAE, self).__init__()
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, latent_dim * 2)  # mean and log-variance
        )
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, input_dim),
            nn.Sigmoid()
        )

    def reparameterize(self, mu, logvar):
        std = torch.exp(0.5 * logvar)
        eps = torch.randn_like(std)
        return mu + eps * std

    def forward(self, x):
        h = self.encoder(x)
        mu, logvar = h.chunk(2, dim=-1)
        z = self.reparameterize(mu, logvar)
        return self.decoder(z), mu, logvar

def loss_function(recon_x, x, mu, logvar):
    BCE = nn.functional.binary_cross_entropy(recon_x, x, reduction='sum')
    KLD = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())
    return BCE + KLD

# データの生成
torch.manual_seed(0)
data = torch.randn(1000, 20)

# モデルの初期化とトレーニング
vae = VAE(input_dim=20, hidden_dim=50, latent_dim=10)
optimizer = optim.Adam(vae.parameters(), lr=1e-3)
num_epochs = 50

for epoch in range(num_epochs):
    vae.train()
    optimizer.zero_grad()
    recon_data, mu, logvar = vae(data)
    loss = loss_function(recon_data, data, mu, logvar)
    loss.backward()
    optimizer.step()

    if (epoch + 1) % 10 == 0:
        print(f'Epoch {epoch+1}, Loss: {loss.item():.4f}')

print("Training complete!")

# 生成結果の表示
vae.eval()
with torch.no_grad():
    sample = torch.randn(64, 10)
    generated_data = vae.decoder(sample).numpy()

plt.figure(figsize=(8, 6))
plt.imshow(generated_data, aspect='auto', cmap='viridis')
plt.colorbar()
plt.title("Generated Data (VAE)")
plt.xlabel("Features")
plt.ylabel("Samples")
plt.show()

# 必要なライブラリのインポート
import numpy as np

# 4.1.1 正方行列の生成
print("4.1.1 Square Matrices")
A = np.array([[1, 2], [3, 4]])
print("Square Matrix A:\n", A)

# 4.1.2 対称行列の生成
print("\n4.1.2 Symmetric Matrices")
B = np.array([[1, 2], [2, 1]])
print("Symmetric Matrix B:\n", B)

# 4.1.3 半正定値行列と正定値行列の生成
print("\n4.1.3 Positive Semidefinite and Positive Definite Matrices")
C = np.array([[2, -1], [-1, 2]])
D = np.array([[2, 1], [1, 2]])
print("Positive Semidefinite Matrix C:\n", C)
print("Positive Definite Matrix D:\n", D)

# 半正定値行列の検証
def is_positive_semidefinite(matrix):
    eigenvalues = np.linalg.eigvals(matrix)
    return np.all(eigenvalues >= 0)

# 正定値行列の検証
def is_positive_definite(matrix):
    eigenvalues = np.linalg.eigvals(matrix)
    return np.all(eigenvalues > 0)

print("Is C positive semidefinite?", is_positive_semidefinite(C))
print("Is D positive definite?", is_positive_definite(D))

# 4.2 行列全体からなるユークリッド空間
print("\n4.2 Euclidean Space Comprising Entire Matrices")

# 4.2.1 内積とノルムの計算
print("\n4.2.1 Inner Product and Norm")
E = np.array([[1, 2], [2, 1]])
F = np.array([[3, 4], [4, 3]])
inner_product = np.trace(np.dot(E.T, F))
norm_E = np.linalg.norm(E, 'fro')
norm_F = np.linalg.norm(F, 'fro')
print("Inner Product of E and F:", inner_product)
print("Frobenius Norm of E:", norm_E)
print("Frobenius Norm of F:", norm_F)

# 4.2.2 行列点列の収束性の検証
print("\n4.2.2 Convergence of Matrix Sequences")
def matrix_sequence_convergence(sequence, tol=1e-5):
    for i in range(1, len(sequence)):
        if np.linalg.norm(sequence[i] - sequence[i-1], 'fro') > tol:
            return False
    return True

# 行列点列の生成と収束性の検証
sequence = [np.eye(2) / (i + 1) for i in range(10)]
print("Matrix sequence converges:", matrix_sequence_convergence(sequence))

# 例としての機械学習への応用
# データの生成と正定値カーネル行列の計算 (RBFカーネル)
from sklearn.metrics.pairwise import rbf_kernel

X = np.random.randn(100, 2)
gamma = 0.1
K = rbf_kernel(X, gamma=gamma)

print("\nGenerated RBF Kernel Matrix K is positive definite:", is_positive_definite(K))

# 必要なライブラリのインポート
import numpy as np
import scipy.linalg as la

# 6.1 固有値問題 (Eigenvalue Problem)
print("6.1 Eigenvalue Problem")
A = np.array([[4, -2], [1, 1]])
eigenvalues, eigenvectors = np.linalg.eig(A)
print("Matrix A:\n", A)
print("Eigenvalues of A:", eigenvalues)
print("Eigenvectors of A:\n", eigenvectors)

# 6.2 エルミート行列の固有値問題 (Eigenvalue Problem of Hermitian Matrices)
print("\n6.2 Eigenvalue Problem of Hermitian Matrices")
# エルミート行列の生成
H = np.array([[2, 1+1j], [1-1j, 3]])
eigenvalues_H, eigenvectors_H = np.linalg.eigh(H)  # Hermitianの場合、eighを使用
print("Hermitian Matrix H:\n", H)
print("Eigenvalues of H:", eigenvalues_H)
print("Eigenvectors of H:\n", eigenvectors_H)

# 6.3 一般化固有値問題 (Generalized Eigenvalue Problem)
print("\n6.3 Generalized Eigenvalue Problem")
# 一般化固有値問題の行列を生成
A = np.array([[1, 2], [3, 4]])
B = np.array([[2, 0], [0, 1]])
eigenvalues_generalized, eigenvectors_generalized = la.eig(A, B)
print("Matrix A:\n", A)
print("Matrix B:\n", B)
print("Generalized Eigenvalues of A and B:", eigenvalues_generalized)
print("Generalized Eigenvectors of A and B:\n", eigenvectors_generalized)

# 固有値分解の可視化
import matplotlib.pyplot as plt

# 固有値のプロット
plt.figure(figsize=(12, 6))

# 固有値問題
plt.subplot(1, 3, 1)
plt.scatter(np.real(eigenvalues), np.imag(eigenvalues), color='red')
plt.title("Eigenvalues of A")
plt.xlabel("Real part")
plt.ylabel("Imaginary part")
plt.grid(True)

# エルミート行列の固有値問題
plt.subplot(1, 3, 2)
plt.scatter(np.real(eigenvalues_H), np.imag(eigenvalues_H), color='blue')
plt.title("Eigenvalues of Hermitian Matrix H")
plt.xlabel("Real part")
plt.ylabel("Imaginary part")
plt.grid(True)

# 一般化固有値問題
plt.subplot(1, 3, 3)
plt.scatter(np.real(eigenvalues_generalized), np.imag(eigenvalues_generalized), color='green')
plt.title("Generalized Eigenvalues of A and B")
plt.xlabel("Real part")
plt.ylabel("Imaginary part")
plt.grid(True)

plt.tight_layout()
plt.show()

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt

# Tepliz行列の生成関数
def generate_tepliz_matrix(c, r):
    return la.toeplitz(c, r)

# 例としてのTepliz行列の生成
c = np.array([1, 2, 3, 4])
r = np.array([1, -1, -2, -3])
tepliz_matrix = generate_tepliz_matrix(c, r)

print("Tepliz Matrix:\n", tepliz_matrix)

# 工学応用: 線形システムの解法
# Ax = b を解く
A = tepliz_matrix
b = np.array([1, 2, 3, 4])

# 解法1: 直接解く
x_direct = np.linalg.solve(A, b)
print("\nSolution using np.linalg.solve:\n", x_direct)

# 解法2: ルーチン法 (SciPyを使用)
x_routine = la.solve_toeplitz((c, r), b)
print("\nSolution using la.solve_toeplitz:\n", x_routine)

# 結果のプロット
plt.figure(figsize=(10, 6))
plt.plot(x_direct, 'o-', label='Direct Solution')
plt.plot(x_routine, 's-', label='Routine Solution')
plt.xlabel('Index')
plt.ylabel('Value')
plt.title('Solution of Linear System Using Tepliz Matrix')
plt.legend()
plt.grid(True)
plt.show()