python进行矩阵运算

news2024/12/19 20:14:49

python进行矩阵运算

一、矩阵相乘
二、矩阵求逆
三、矩阵特征值与特征向量求解
四、矩阵约当标准型与转换矩阵求解
五、矩阵奇异值分解
六、矩阵方程组求解

一、矩阵相乘

python中矩阵相乘可以使用numpy实现，也可以使用sympy实现，以numpy实现为例，代码如下：

import numpy as np

a01=np.random.randint(10,size=(4,3))
print(a01)
array([[6, 4, 9],
       [1, 4, 8],
       [8, 0, 0],
       [5, 6, 8]])
a02=np.random.randint(10,size=(3,4))
print(a02)
array([[5, 7, 7, 9],
       [6, 4, 6, 2],
       [5, 2, 8, 2]])

a03=np.dot(a01,a02)
print(a03)
array([[ 99,  76, 138,  80],
       [ 69,  39,  95,  33],
       [ 40,  56,  56,  72],
       [101,  75, 135,  73]])

也可将乘积转换为分数形式，示例代码如下：

import numpy as np
from fractions import Fraction

a01=np.random.randint(10,size=(4,3)).astype(float)
a01[0,0]=0.5
print(a01)
array([[0.5, 4. , 3. ],
       [0. , 4. , 4. ],
       [4. , 3. , 6. ],
       [0. , 2. , 6. ]])
a02=np.random.randint(10,size=(3,4))
print(a02)
array([[5, 7, 7, 9],
       [6, 4, 6, 2],
       [5, 2, 8, 2]])

a03=np.dot(a01,a02)
print(a03)#a03为小数形式
array([[41.5, 25.5, 51.5, 18.5],
       [44. , 24. , 56. , 16. ],
       [68. , 52. , 94. , 54. ],
       [42. , 20. , 60. , 16. ]])

a04=np.vectorize(lambda x: Fraction.from_float(x).limit_denominator())(a03)
print(a04)#a04为分数形式
array([[Fraction(83, 2), Fraction(51, 2), Fraction(103, 2), Fraction(37, 2)],
       [Fraction(44, 1), Fraction(24, 1), Fraction(56, 1), Fraction(16, 1)],
       [Fraction(68, 1), Fraction(52, 1), Fraction(94, 1), Fraction(54, 1)],
       [Fraction(42, 1), Fraction(20, 1), Fraction(60, 1), Fraction(16, 1)]], dtype=object)

二、矩阵求逆

分别使用numpy和sympy实现，代码如下：

import numpy as np
from fractions import Fraction

#numpy实现
a01=np.random.randint(5,size=(2,2))
print(a01)
array([[4, 0],
       [3, 4]])
a02 = np.linalg.inv(a01.astype(float))#求逆矩阵
print(a02)
array([[ 0.25  ,  0.    ],
       [-0.1875,  0.25  ]])

a03=np.vectorize(lambda x: Fraction.from_float(x).limit_denominator())(a02)
print(a03)#分数形式
array([[Fraction(1, 4), Fraction(0, 1)],
       [Fraction(-3, 16), Fraction(1, 4)]], dtype=object)

#sympy实现
a04=Matrix([[4,0],[3,4]])
print(a04)
Matrix([
[4, 0],
[3, 4]])
print(a04.inv())
Matrix([
[  1/4,   0],
[-3/16, 1/4]])

从求解结果分析，sympy会自动给出结果的简洁形式，如果是分数则用分数表示，在展示效果上优于numpy。

三、矩阵特征值与特征向量求解

分别使用numpy和sympy实现，代码如下：

import numpy as np
from sympy import Matrix

#numpy实现
a01=np.array([[-1,2,0],[0,3,0],[2,1,-1]])
print(a01)
array([[-1,  2,  0],
       [ 0,  3,  0],
       [ 2,  1, -1]])
eigenvalue, featurevector = np.linalg.eig(a01)
print(eigenvalue)
array([-1., -1.,  3.])
print(featurevector)
array([[ 0.00000000e+00,  1.11022302e-16,  4.08248290e-01],
       [ 0.00000000e+00,  0.00000000e+00,  8.16496581e-01],
       [ 1.00000000e+00, -1.00000000e+00,  4.08248290e-01]])

#sympy实现
a02=Matrix([[-1,2,0],[0,3,0],[2,1,-1]])
print(a02)
Matrix([
[-1, 2,  0],
[ 0, 3,  0],
[ 2, 1, -1]])
print(a02.eigenvals())
{3: 1, -1: 2}#3: 1意为特征值为3，个数为1个，-1: 2意为特征值为-1，个数为2个
print(a02.eigenvects())
[(-1, 2, [Matrix([
[0],
[0],
[1]])]), (3, 1, [Matrix([
[1],
[2],
[1]])])]

从求解结果看，sympy的求解结果比numpy求解结果更为直观。

四、矩阵约当标准型与转换矩阵求解

矩阵约当标准型与转换矩阵使用numpy求解较为繁琐，因此使用sympy求解，代码如下：

from sympy import Matrix

a01=Matrix([[-1,2,0],[0,3,0],[2,1,-1]])
print(a01)
Matrix([
[-1, 2,  0],
[ 0, 3,  0],
[ 2, 1, -1]])
p_mat,j_mat=a01.jordan_form()#p_mat为转换矩阵，j_mat为约当标准型
print(p_mat)
Matrix([
[0, 1, 1],
[0, 0, 2],
[2, 0, 1]])
print(j_mat)
Matrix([
[-1,  1, 0],
[ 0, -1, 0],
[ 0,  0, 3]])

注意sympy求解的约当标准型的1放在上三角，有些数学教材中将1放在下三角。

五、矩阵奇异值分解

sympy不直接支持矩阵奇异值的求解，因此使用numpy求解，代码如下：

import numpy as np

a02=np.array([[4,0],[3,0],[0,0]])
print(a02)
array([[4, 0],
       [3, 0],
       [0, 0]])
U, S, Vt = np.linalg.svd(a02)
print(U)#左奇异向量组成的酉矩阵
array([[-0.8, -0.6,  0. ],
       [-0.6,  0.8,  0. ],
       [ 0. ,  0. ,  1. ]])
print(S)#奇异值
array([5., 0.])
print(Vt)#右奇异向量组成的酉矩阵
array([[-1., -0.],
       [ 0.,  1.]])

奇异向量组成的酉矩阵并不唯一。

六、矩阵方程组求解

使用sympy求解，代码如下：

from sympy import Matrix, symbols, linsolve

A = Matrix([[2,3,1],[4,2,3], [7,1,-1]])
print(A)
Matrix([
[2, 3,  1],
[4, 2,  3],
[7, 1, -1]])
B = Matrix([[4],[17],[1]])
print(B)
Matrix([
[ 4],
[17],
[ 1]])
x, y, z = symbols('x y z')
result=linsolve((A,B),x,y,z)
print(result)
{(1, -1, 5)}

A = Matrix([[1,1,1,1],[4,3,5,-1], [2,1,3,-3]])
print(A)
Matrix([
[1, 1, 1,  1],
[4, 3, 5, -1],
[2, 1, 3, -3]])
B = Matrix([[-1],[-1],[1]])
print(B)
Matrix([
[-1],
[-1],
[ 1]])
x, y, z = symbols('x y z')
result=linsolve((A,B),x,y,z)
print(result)
{(2 - 2*z, z - 3, z)}