python-numpy(3)-线性代数

一、方程求解

参考资料
对于Ax = b 这种方程：

1. np.linalg.inv(A).dot(B)
2. np.linalg.solve(A,b)

1.1 求解多元一次方程一个直观的例子

# AX=B
# X = A^(-1)*B 
A = np.array([
    [7, 3, 0, 1], 
    [0, 1, 0, -1], 
    [1, 0, 6, -3], 
    [1, 1, -1, -1]])
B = np.array([
    8, 
    6, 
    -3, 
    1])
X = np.linalg.inv(A).dot(B)

print("x是:{} ".format(X[0]))
print("y是:{} ".format(X[1]))
print("z是:{} ".format(X[2]))
print("k是:{} ".format(X[3]))

1.2. 方程求解具体步骤

1.2.1 创建矩阵

A = np.mat("0 1 2;1 0 3;4 -3 8")
print (A)
#[[ 0 1 2]
# [ 1 0 3]
# [ 4 -3 8]]

1.2.2 求解逆矩阵

# 使用inv函数计算逆矩阵
inv = np.linalg.inv(A)
print (inv)
#[[-4.5 7. -1.5]
# [-2. 4. -1. ]
# [ 1.5 -2. 0.5]]

验证是否是逆矩阵

print (A * inv)
#[[ 1. 0. 0.]
# [ 0. 1. 0.]
# [ 0. 0. 1.]]

1.2.3 完整的求解过程

# Bx = b 
#创建矩阵和数组
B = np.mat("1 -2 1;0 2 -8;-4 5 9")
b = np.array([0,8,-9])

# 调用solve函数求解线性方程
x = np.linalg.solve(B,b)
print (x)
#[ 29. 16. 3.]

# 使用dot函数检查求得的解是否正确
print (np.dot(B , x))
# [[ 0. 8. -9.]]

二、常用的定义

2.1 特征值和特征向量

• 特征值（eigenvalue）即方程 Ax = ax 的根，是一个标量。其中，A 是一个二维矩阵，x 是一个一维向量。特征向量（eigenvector）是关于特征值的向量
• numpy.linalg模块中，eigvals函数可以计算矩阵的特征值，而eig函数可以返回一个包含特征值和对应的特征向量的元组

# 创建一个矩阵
C = np.mat("3 -2;1 0")

# 调用eigvals函数求解特征值
c0 = np.linalg.eigvals(C)
print (c0)
# [ 2. 1.]

# 使用eig函数求解特征值和特征向量 (该函数将返回一个元组，按列排放着特征值和对应的特征向量，其中第一列为特征值，第二列为特征向量)
c1,c2 = np.linalg.eig(C)
print (c1)
# [ 2. 1.]
print (c2)
#[[ 0.89442719 0.70710678]
# [ 0.4472136 0.70710678]]

# 使用dot函数验证求得的解是否正确
for i in range(len(c1)):
    print ("left:",np.dot(C,c2[:,i]))
    print ("right:",c1[i] * c2[:,i])
#left: [[ 1.78885438]
# [ 0.89442719]]
#right: [[ 1.78885438]
# [ 0.89442719]]
#left: [[ 0.70710678]
# [ 0.70710678]]
#right: [[ 0.70710678]
# [ 0.70710678]]

使用数组函数创建数组

a = np.array([[1, -2j], [2j, 5]])

print("Array is :",a)
# 使用 with() 函数计算特征值
c, d = np.linalg.eigh(a)

print("Eigen value is :", c)
print("Eigen value is :", d)

2.2 奇异值分解

参考资料

• SVD（Singular Value Decomposition，奇异值分解）是一种因子分解运算，将一个矩阵分解为3个矩阵的乘积
• numpy.linalg模块中的svd函数可以对矩阵进行奇异值分解。该函数返回3个矩阵——U、Sigma和V，其中U和V是正交矩阵，Sigma包含输入矩阵的奇异值。

# 分解矩阵
D = np.mat("4 11 14;8 7 -2")
# 使用svd函数分解矩阵
U,Sigma,V = np.linalg.svd(D, full_matrices=False)
print ("U:",U)
#U: [[-0.9486833 -0.31622777]
# [-0.31622777 0.9486833 ]]
print ("Sigma:",Sigma)
#Sigma: [ 18.97366596 9.48683298]
print ("V",V)
#V [[-0.33333333 -0.66666667 -0.66666667]
# [ 0.66666667 0.33333333 -0.66666667]]
# 结果包含等式中左右两端的两个正交矩阵U和V，以及中间的奇异值矩阵Sigma

# 使用diag函数生成完整的奇异值矩阵。将分解出的3个矩阵相乘
print (U * np.diag(Sigma) * V)
#[[ 4. 11. 14.]
# [ 8. 7. -2.]]

2.3 逆矩阵

• 使用numpy.linalg模块中的pinv函数进行求解,
• 注：inv函数只接受方阵作为输入矩阵，而pinv函数则没有这个限制

2.3.1 普通逆矩阵

# 使用inv函数计算逆矩阵
inv = np.linalg.inv(A)
print (inv)
#[[-4.5 7. -1.5]
# [-2. 4. -1. ]
# [ 1.5 -2. 0.5]]

2.3.2 广义逆矩阵

# 创建一个矩阵
E = np.mat("4 11 14;8 7 -2")
# 使用pinv函数计算广义逆矩阵
pseudoinv = np.linalg.pinv(E)
print (pseudoinv)
#[[-0.00555556 0.07222222]
# [ 0.02222222 0.04444444]
# [ 0.05555556 -0.05555556]]

# 将原矩阵和得到的广义逆矩阵相乘
print (E * pseudoinv)
#[[ 1.00000000e+00 -5.55111512e-16]
# [ 0.00000000e+00 1.00000000e+00]]

2.4 行列式

• numpy.linalg模块中的det函数可以计算矩阵的行列式

# 计算矩阵的行列式
F = np.mat("3 4;5 6")
# 使用det函数计算行列式
print (np.linalg.det(F))
# -2.0

2.5 范数

• 顾名思义，linalg=linear+algebra，normnorm则表示范数，首先需要注意的是范数是对向量（或者矩阵）的度量，是一个标量（scalar）：
• 首先help(np.linalg.norm)查看其文档：norm(x, ord=None, axis=None, keepdims=False)
  *

a=[[1,2,0],[-1,2,-1],[0,1,1]]#写入矩阵
A=np.array(a)#转换成np.array格式
a1=np.linalg.norm(A, ord=1)
a2=np.linalg.norm(A, ord=2)
a2_2=np.sum(np.abs(A) ** 2) ** 0.5 

a3=np.linalg.norm(A, ord=np.Inf)
a4=np.linalg.norm(A, ord=-np.Inf)
print("第一范数为",a1)
print("第二范数为",round(a2, 2))
print("第二范数_自定义函数为",round(a2_2, 2))
print("无穷大范数为",a3)
print("无穷小范数a4为",a3)

#法二：第二范数：先求最大特征值再求
 
#先得到A的转置矩阵与A相乘的矩阵d
b=np.transpose(A)
B=np.array(b)
d=np.matmul(B,A)
 
#根据特征多项式得一元三次方程求解
import sympy as sp # 导入sympy包
x=sp.Symbol('x')
f=x**3-13*(x**2)+38*x-25   #见前面运算过程
x=sp.solve(f)  #得到的解析解比较复杂，故后续转换成浮点数
 
#取四位小数输出
for i in range(0,3):
    x[i]=round(x[i].evalf(),4) #求出表达式的浮点数
    print(x[i])
maxX=x[0]
for i in x:   ## 求最大值
    if i > maxX:
        maxX= i
print("第二范数为",round(maxX**0.5,4)) 

2.6 矩阵的秩

import numpy as np
A = np.array([
        [1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
        [0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0],
        [5, 1, 2, 4, 0, 0, 1, 0, 0, 0, 0, 0],
        [5, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
        [1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
        [0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
        [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0],
        [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1]
])
B = np.array(
    [
        [1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 10],
        [0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 15],
        [5, 1, 2, 4, 0, 0, 1, 0, 0, 0, 0, 0, 100],
        [5, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 60],
        [1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,1],
        [0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,4],
        [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0,2],
        [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1,10]
]
)
np.linalg.matrix_rank(A)#返回矩阵的秩
np.linalg.matrix_rank(B)#返回矩阵的秩

2.7 矩阵的迹

方阵的迹就是主对角元素之和：

trace(a, offset=0, axis1=0, axis2=1, dtype=None, out=None):

import numpy as np

x = np.array([[1, 2, 3], [3, 4, 5], [6, 7, 8]])
print(x)
# [[1 2 3]
#  [3 4 5]
#  [6 7 8]]

y = np.array([[5, 4, 2], [1, 7, 9], [0, 4, 5]])
print(y)
# [[5 4 2]
#  [1 7 9]
#  [0 4 5]]

print(np.trace(x))  # A的迹等于A.T的迹
# 13
print(np.trace(np.transpose(x)))
# 13

print(np.trace(x + y))  # 和的迹 等于 迹的和
# 30
print(np.trace(x) + np.trace(y))
# 30

2.8 两个数组的矩阵积

numpy.matmul 函数返回两个数组的矩阵乘积。

• 虽然它返回二维数组的正常乘积，
• 如果任一参数的维数大于2，则将其视为存在于最后两个索引的矩阵的栈，并进行相应广播。
• 另一方面，如果任一参数是一维数组，则通过在其维度上附加 1 来将其提升为矩阵，并在乘法之后被去除。

import numpy as np 
 
a = [[1,0],[0,1]] 
b = [[4,1],[2,2]] 
print (np.matmul(a,b))

a = [[1,0],[0,1]] 
b = [1,2] 
print (np.matmul(a,b))
print (np.matmul(b,a))

a = np.arange(8).reshape(2,2,2) 
b = np.arange(4).reshape(2,2) 
print (np.matmul(a,b))

2.9 两个向量的内积

numpy.inner() 函数返回一维数组的向量内积。对于更高的维度，它返回最后一个轴上的和的乘积。

print (np.inner(np.array([1,2,3]),np.array([0,1,0])))
# 等价于 1*0+2*1+3*0

a = np.array([[1,2], [3,4]]) 
b = np.array([[11, 12], [13, 14]])  
"""
1*11+2*12, 1*13+2*14 
3*11+4*12, 3*13+4*14
"""
print ('内积：')
print (np.inner(a,b))