用python进行数值分析 - 矩阵基础

import numpy as np
from scipy import linalg

# 求解多元线性方程组必定涉及到矩阵，也必定涉及到矩阵的一些常用概念和操作
# 比如矩阵的特征值、特征向量，求逆矩阵，矩阵分解等
# 以下是常用的几个矩阵的操作
# 参考资料: https://docs.scipy.org/doc/scipy/tutorial/linalg.html
# 矩阵的特征值和特征向量
# 定义:
# A 为nxn的矩阵，x为非零向量，若存在数λ使Ax = λx有非平凡解x, 则称λ为A的特征值, x称为对应于λ的特征向量。
# 非平凡解就是不全为零的一个向量
print("===============求矩阵的特征值和特征向量===================")
A = np.array([[1, 5, 2], [2, 4, 1], [3, 6, 2]])
print('数组A:\n', A)
λ, v = linalg.eig(A)
print('数组A的特征值是: \n', λ)
print('数组A的特征向量是: \n', v)
# 判定所求的结果是否准确呢？
part1 = A.dot(v)
part2 = λ*v
print('等号左边的矩阵为: \n', part1)
print('等号右边的矩阵为: \n', part2)

print("================矩阵的逆(numpy)==================")
# 使用numpy求矩阵的逆
a = np.mat('[1 2; 3 4]')
print(a.I)
b = np.mat('[-2 1; 1.5 -0.5]')
print(a*b)


print("================矩阵的逆(scipy)==================")
# 使用scipy求矩阵的逆
# inv API doc: 
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.inv.html
sa = np.array([[1., 2.], [3., 4.]])
sai = linalg.inv(sa, check_finite=True)
print(sai)
print(np.dot(sa, sai)) #并不是期望的结果！！
# 手动输入结果后的结果验证
bi = np.array([[-2., 1.], [1.5, -0.5]])
si = np.dot(sa, bi)
print(si)


print("================矩阵LU分解==================")
# LU分解一般主要针对方阵进行
# 初等行变换即可进行，L表示下三角矩阵，U表示上三角矩阵
# 一般的分解条件是: 可逆的方阵
A = np.array([[1, 2, -1], [2, 1, -2], [-3, 1, 1]])
print("矩阵A: \n", A)
Ai = linalg.inv(A)
print("矩阵A的逆矩阵为: \n", Ai)
# lu API 文档: 
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu.html
l, u = linalg.lu(A, permute_l=True) # 如果期望结果只包含L和U，参数permute_l要设置成True
print("矩阵A的LU分解 - l:\n", l)
print("矩阵A的LU分解 - u:\n", u)
A1 = l@u
print("验证PLU分解的结果: \n", A1)
p, l, u = linalg.lu(A)
print("矩阵A的LU分解 - p:\n", p)
print("矩阵A的LU分解 - l:\n", l)
print("矩阵A的LU分解 - u:\n", u)
A2 = p@l@u
print("验证PLU分解的结果: \n", A2)


print("================矩阵楚列斯基分解==================")
# 楚列斯基分解定理：
# 如果A是nxn对称正定矩阵，则存在上三角nxn矩阵R满足A = Rt*R (Rt表示R的转置)
# API 文档: 
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.cholesky.html
M = np.array([[4, -2, 2], [-2, 2, -4], [2, -4, 11]])
print("矩阵M为:\n", M)
L = linalg.cholesky(M, lower=True)
print("矩阵M的楚列斯基分解得到的下三角矩阵为:\n", L)
H = linalg.cholesky(M, lower=False)
print("矩阵M的楚列斯基分解得到的上三角矩阵为:\n", H)
# 对称矩阵很容易判定，即Mt = M, 怎么判定一个矩阵是不是正定矩阵呢？
# 如果nxn矩阵A是对称矩阵，则当且仅当A的所有特征值是正数时A为正定矩阵
# eigvals API 文档: https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigvals.html
λ1 = linalg.eigvals(M)
print('数组M的特征值是: \n', λ1)
for n in λ1:
    print('特征值大于0吗?: ', n>0) 
M1 = H.T@H
print('矩阵M的楚列斯基分解的上三角验证: \n', M1)
M2 = L@L.T
print('矩阵M的楚列斯基分解的下三角验证: \n', M2)