书的作者
文章目录
- 算术乘除：先乘除，后加减，括号内先算
- 基本的乘法运算
- 计算阶乘
- 基本除法
- 向量的乘法：标量乘法，向量内积，逐项积
- 标量乘法
- 向量的内积
- 对于inner和dot的实现方式的探究
- 逐项积
- dot的计算过程
- 逐项积和向量内积的区别
- 矩阵乘法：最重要的线性代数运算规则
- 矩阵乘法规则
- 矩阵乘法的形态
- 向量相乘
- 矩阵乘法的第一个视角
- Streamlit程序
- 矩阵乘法的第二视角
- 代码：矩阵乘法的热图和符号的可视化
- 矩阵除法：计算逆矩阵
- 广义逆矩阵
- 代码inv求非奇异方阵的逆
使用到的函数
重点：
向量乘法的三个大类
矩阵乘法的基础
算术乘除：先乘除，后加减，括号内先算
基本的乘法运算
num1 = 2
num2 = 3
# add two numbers
prod = num1*num2
# display the computation
print('The product of {0} and {1} is {2}'.format(num1, num2, prod))
计算阶乘
判断限制条件，循环累乘
num = int(input("Enter an integer: "))
factorial = 1
# check if the number is negative, positive or zero
if num < 0:
print("Factorial does not exist for negative numbers")
elif num == 0:
print("The factorial of 0 is ", factorial)
else:
for i in range(1,num + 1):
factorial = factorial*i
print("The factorial of",num," is ",factorial)
直接调用scipy的阶乘算法
import scipy
scipy.special.factorial(5)
利用 numpy.linspace(1, 10, 10) 产生 1 ~ 10 这 10 个自然数，然后利用 numpy.cumprod() 函数来求累计乘积。
import numpy as np
a_i = np.linspace(1,10,10)
print(a_i)
#累加是这个np.cumsum(a_i),上一章出现过
#a_i = np.linspace(1,10,10)
#累乘
a_i_cumprod = np.cumprod(a_i)
# 实现比较简单，使用阶乘，把每次的值记录下来就行
#suppress可以让数据展现完整
np.set_printoptions(suppress=True)
print(a_i_cumprod)
#使用函数实现的累乘
import numpy as np
a_i = np.linspace(1,10,10)
print(a_i)
#累乘
a_i_cumprod = np.cumprod(a_i)
# 实现比较简单，使用阶乘，把每次的值记录下来就行
#suppress可以让数据展现完整
np.set_printoptions(suppress=True)
print(a_i_cumprod)
[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]
[ 1. 2. 6. 24. 120. 720. 5040. 40320. 362880. 3628800.]
基本除法
除数
num1 = 6
num2 = 3
# add two numbers
division = num1/num2
# display the computation
print('The division of {0} over {1} is {2}'.format(num1, num2, division))
余数
num1 = 7
num2 = 2
# add two numbers
remainder = num1%num2
# display the computation
print('The remainder of {0} over {1} is {2}'.format(num1, num2, remainder))
倒数
def reciprocal(x):
return 1 / x
#幂运算
def reciprocal(x):
return x ** -1
#对于需要精确分数表示的情况，使用fraction模块
from fractions import Fraction
def reciprocal(x):
return Fraction(1, x)
向量的乘法：标量乘法，向量内积，逐项积
三种向量乘法
（1）标量乘法(scalar multiplication)
（2）向量内积(inner product)
（3）逐项积(piecewise product)
标量乘法
标量乘法运算中，标量乘向量结果还是向量，相当于缩放。
import numpy as np
# define a column vector
a_col = np.array([[1], [2], [3]])
b_col = 2*a_col
# define a row vector
a_row = np.array([[1, 2, 3]])
b_row = 2*a_row
# define a matrix
A = np.array([[1, 2, 3],
[4, 5, 6]])
B = 2*A
向量的内积
定律
import numpy as np
# row vector dot product
a_row = np.array([[1, 2, 3]])
b_row = np.array([[4, 3, 2]])
a_dot_b = np.inner(a_row, b_row)
print(a_dot_b)
print(np.inner(a_row[:], b_row[:]))
print(np.sum(a_row * b_row))
#%% column vector dot product
a_col = np.array([[1], [2], [3]])
b_col = np.array([[-1], [0], [1]])
a_dot_b = np.inner(a_col, b_col)
print(a_dot_b) # tensor product
print(np.sum(a_col * b_col))
import numpy as np
#[1,3]
a_row = np.array([[1, 2, 3]])
b_row = np.array([[4, 3, 2]])
#向量内积
a_dot_b = np.inner(a_row, b_row)
print(a_dot_b)
#[[16]]
print(np.inner(a_row[:], b_row[:]))
#[[16]]
#先逐项积再求和
print(np.sum(a_row * b_row))
#16
#%% column vector dot product
#[3,1]
a_col = np.array([[1], [2], [3]])
b_col = np.array([[-1], [0], [1]])
a_dot_b = np.inner(a_col, b_col)
#变成了广播机制类似的
print(a_dot_b) # tensor product
#[[-1 0 1]
#[-2 0 2]
#[-3 0 3]]
print(np.sum(a_col * b_col))
#2
-
np.inner函数的行为：
np.inner(a, b)计算两个数组的内积。它的行为会根据输入数组的维度而变化。 -
对于一维数组：
如果a_col和b_col都是一维数组，np.inner会计算它们的点积，不需要额外的对齐。 -
对于多维数组：
- 如果
a_col和b_col是多维数组，np.inner会对最后一个轴进行求和。 - 它不会自动对齐其他维度。两个数组除了最后一个轴外，其他轴的大小必须相同。
- 如果
-
列向量的情况：
- 如果
a_col和b_col是列向量（shape 为 (n, 1) 的二维数组），np.inner不会自动将它们视为一维数组。 - 在这种情况下，您可能需要先将它们转换为一维数组（通过
.flatten()或.ravel()），或者使用np.dot()函数。
- 如果
上面
a_col = np.array([[1], [2], [3]])
b_col = np.array([[-1], [0], [1]])
就是列向量，无法正确计算结果，
a_col.flatten()
#array([1, 2, 3]) ,shape(3,)
np.inner 的行为：
- 对于多维数组，np.inner 对最后一个轴进行求和，而对其他轴进行广播。
- 在这个案例中，a_col 和 b_col 都是形状为 (3, 1) 的二维数组。
实际发生的计算：
- np.inner 会将这两个数组视为一组一维数组。
- 它会对每对一维数组（在这里，每个一维数组只有一个元素）计算内积。
- 结果是一个 3x3 的矩阵，而不是我们可能期望的单个数字。
示例：
import numpy as np
# 一维数组 - 正常工作
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
print(np.inner(a, b)) # 输出: 32
# 二维列向量 - 可能不是您想要的结果
a_col = np.array([[1], [2], [3]])
b_col = np.array([[4], [5], [6]])
print(np.inner(a_col, b_col)) # 输出: [[32]]
# 正确处理列向量的方法
print(np.inner(a_col.flatten(), b_col.flatten())) # 输出: 32
# 或者
print(np.dot(a_col.T, b_col).item()) # 输出: 32
总结：np.inner 不会自动对齐所有维度。对于列向量，您可能需要手动将它们转换为一维数组或使用其他函数（如 np.dot）来获得预期的结果。
对于inner和dot的实现方式的探究
矩阵和向量的乘法是一种常见的矩阵运算，也被称为矩阵-向量乘积。它的定义是，将一个 m × n m \times n m×n 的矩阵 A A A 乘以一个 n × 1 n \times 1 n×1 的列向量 x \mathbf{x} x，得到一个 m × 1 m \times 1 m×1 的列向量 y \mathbf{y} y，其中 m m m 表示矩阵 A A A 的行数， n n n 表示矩阵 A A A 的列数。
具体来说，矩阵和向量的乘法可以使用矩阵的行向量和向量的列向量逐一相乘的方式来实现。可以将矩阵 A A A 拆分为它的行向量 A i , : A_{i,:} Ai,:，将向量 x \mathbf{x} x 拆分为它的列向量 x : , 1 \mathbf{x}_{:,1} x:,1，然后对它们逐一相乘并相加，得到矩阵和向量的乘积 y \mathbf{y} y 的每个元素 y i y_i yi：
y
i
=
∑
j
=
1
n
A
i
,
j
x
j
,
i
=
1
,
…
,
m
y_i = \sum_{j=1}^{n} A_{i,j} x_j ,\quad i=1,\ldots,m
yi=j=1∑nAi,jxj,i=1,…,m
ndim返回数组的维度
逐项积
import numpy as np
a = np.array([[1, 2, 3]])
b = np.array([[4, 5, 6]])
# calculate element-wise product of row vectors
a_times_b = a*b
A = np.array([[1, 2, 3],
[4, 5, 6]])
B = np.array([[1, 2, 3],
[-1, 0, 1]])
# calculate element-wise product of matrices
A_times_B = A*B
dot的计算过程
def dot(a, b):
if a.ndim == 1 and b.ndim == 1:
# 一维数组的内积，相当于
result = 0
for i in range(len(a)):
result += a[i] * b[i]
return result
elif a.ndim == 2 and b.ndim == 1:
# 矩阵和向量的乘法
#最后的结果是 a的行阶的向量
result = np.zeros(a.shape[0])
for i in range(a.shape[0]):
#行向量的每个值和b的列向量相乘再相加
for j in range(a.shape[1]):
result[i] += a[i, j] * b[j]
return result
elif a.ndim == 1 and b.ndim == 2:
# 向量和矩阵的乘法，结果为列向量
result = np.zeros(b.shape[1])
for i in range(b.shape[1]):
for j in range(b.shape[0]):
result[i] += a[j] * b[j, i]
return result
elif a.ndim == 2 and b.ndim == 2:
# 矩阵和矩阵的乘法
#a的行，b的列建立一个空数组
result = np.zeros((a.shape[0], b.shape[1]))
for i in range(a.shape[0]):
for j in range(b.shape[1]):
for k in range(a.shape[1]):
result[i, j] += a[i, k] * b[k, j]
return result
逐项积和向量内积的区别
-
逐项积（Element-wise Product）：
- 定义：两个相同维度的向量对应元素相乘。
- 结果：得到一个与原向量相同维度的新向量。
- 数学表示：对于向量 a = [a₁, a₂, …, aₙ] 和 b = [b₁, b₂, …, bₙ]，
逐项积为 [a₁b₁, a₂b₂, …, aₙbₙ] - Python/NumPy 实现：使用 * 运算符或 np.multiply()
import numpy as np a = np.array([1, 2, 3]) b = np.array([4, 5, 6]) element_wise_product = a * b # 或 np.multiply(a, b) print(element_wise_product) # 输出：[4 10 18] -
向量内积（Dot Product 或 Inner Product）：
- 定义：两个向量对应元素相乘后求和。
- 结果：得到一个标量（单个数值）。
- 数学表示：对于向量 a = [a₁, a₂, …, aₙ] 和 b = [b₁, b₂, …, bₙ]，
内积为 a·b = a₁b₁ + a₂b₂ + … + aₙbₙ - Python/NumPy 实现：使用 np.dot() 或 @运算符
import numpy as np a = np.array([1, 2, 3]) b = np.array([4, 5, 6]) dot_product = np.dot(a, b) # 或 a @ b print(dot_product) # 输出：32
主要区别：
-
结果维度：
- 逐项积：结果是一个向量
- 内积：结果是一个标量
-
数学意义：
- 逐项积：保留了每个位置的乘积信息
- 内积：提供了两个向量的相似度或投影信息
-
应用场景：
- 逐项积：常用于元素级的操作，如在神经网络中应用激活函数
- 内积：常用于计算向量间的相似度、投影等，如在机器学习中计算权重和特征的加权和
-
泛化到高维：
- 逐项积：可以直接扩展到多维数组（张量）
- 内积：在矩阵上的扩展是矩阵乘法
矩阵乘法：最重要的线性代数运算规则
矩阵乘法规则
不满足交换律
A
m
×
p
B
p
×
n
≠
B
p
×
n
A
m
×
p
A_{m\times p}B_{p\times n}\neq B_{p\times n}A_{m\times p}
Am×pBp×n=Bp×nAm×p
import numpy as np
A = np.array([[1, 2],
[3, 4]])
B = np.array([[4, 2],
[3, 1]])
# calculate matrix multiplication
C = A@B
矩阵乘法的形态
向量相乘
矩阵乘法的第一个视角
以二维矩阵为例：
A视为行向量的组合，B视为列向量的组合。
Streamlit程序
#运行方式 cmd执行 streamlit run Streamlit_Bk3_Ch2_10.py
###############
# Authored by Weisheng Jiang
# Book 3 | From Basic Arithmetic to Machine Learning
# Published and copyrighted by Tsinghua University Press
# Beijing, China, 2022
###############
import numpy as np
import streamlit as st
import seaborn as sns
from matplotlib import pyplot as plt
#%%
def bmatrix(a):
"""Returns a LaTeX bmatrix
:a: numpy array
:returns: LaTeX bmatrix as a string
"""
if len(a.shape) > 2:
raise ValueError('bmatrix can at most display two dimensions')
lines = str(a).replace('[', '').replace(']', '').splitlines()
rv = [r'\begin{bmatrix}']
rv += [' ' + ' & '.join(l.split()) + r'\\' for l in lines]
rv += [r'\end{bmatrix}']
return '\n'.join(rv)
#%%
with st.sidebar:
st.latex(r'C_{m\times n} = A_{m\times p} B_{p\times n}')
rows_A = st.slider('Number of rows in A:',
min_value = 1,
max_value = 9,
value = 5,
step = 1)
cols_A = st.slider('Number of columns in A:',
min_value = 1,
max_value = 9,
value = 5,
step = 1)
rows_B = st.slider('Number of rows in B:',
min_value = 1,
max_value = 9,
value = 5,
step = 1)
cols_B = st.slider('Number of columns in B:',
min_value = 1,
max_value = 9,
value = 5,
step = 1)
#%% generate A and B using random integer generator
A = np.random.randint(10, size=(rows_A, cols_A))
B = np.random.randint(10, size=(rows_B, cols_B))
st.latex(r'A_{m\times p} = ' + bmatrix(A))
st.latex(r'B_{p\times n} = ' + bmatrix(B))
#%%
try:
C = A@B
st.latex('C = AB = ' + bmatrix(C))
fig, axs = plt.subplots(1, 5, figsize=(12, 3))
plt.sca(axs[0])
ax = sns.heatmap(A,cmap='RdYlBu_r',
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,rows_A+1),
xticklabels=np.arange(1,cols_A+1))
ax.set_aspect("equal")
plt.title('$A$')
plt.yticks(rotation=0)
plt.sca(axs[1])
plt.title('$@$')
plt.axis('off')
plt.sca(axs[2])
ax = sns.heatmap(B,cmap='RdYlBu_r',
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,rows_B+1),
xticklabels=np.arange(1,cols_B+1))
ax.set_aspect("equal")
plt.title('$B$')
plt.yticks(rotation=0)
plt.sca(axs[3])
plt.title('$=$')
plt.axis('off')
plt.sca(axs[4])
ax = sns.heatmap(C,cmap='RdYlBu_r',
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,rows_A+1),
xticklabels=np.arange(1,cols_B+1))
ax.set_aspect("equal")
plt.title('$C$')
plt.yticks(rotation=0)
st.pyplot(fig)
except:
st.write('The number of columns of the first matrix, must equal the number of rows of the second matrix.')
矩阵乘法的第二视角
A视为列向量的组合，是一行，B是一列，是行向量的组合。
按照维度的计算， 应该是一个数的加和。而每个子组合进行计算后的结果不是一个数。而是多个矩阵的叠加。
代码：矩阵乘法的热图和符号的可视化
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
# Repeatability
np.random.seed(7)
# Generate matrix A and B
m = 6
p = 3
n = 4
A = np.random.uniform(-1,1,m*p).reshape(m, p)
B = np.random.uniform(-1,1,p*n).reshape(p, n)
C = A@B
all_max = 1
all_min = -1
#%% matrix multiplication, first perspective
fig, axs = plt.subplots(1, 5, figsize=(12, 3))
plt.sca(axs[0])
ax = sns.heatmap(A,cmap='RdBu_r',vmax = all_max,vmin = all_min,
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,m+1), xticklabels=np.arange(1,p+1))
ax.set_aspect("equal")
plt.title('$A$')
plt.yticks(rotation=0)
plt.sca(axs[1])
plt.title('$@$')
plt.axis('off')
plt.sca(axs[2])
ax = sns.heatmap(B,cmap='RdBu_r',vmax = all_max,vmin = all_min,
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,p+1), xticklabels=np.arange(1,n+1))
ax.set_aspect("equal")
plt.title('$B$')
plt.yticks(rotation=0)
plt.sca(axs[3])
plt.title('$=$')
plt.axis('off')
plt.sca(axs[4])
ax = sns.heatmap(C,cmap='RdBu_r',vmax = all_max,vmin = all_min,
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,m+1), xticklabels=np.arange(1,n+1))
ax.set_aspect("equal")
plt.title('$C$')
plt.yticks(rotation=0)
plt.show()
#%% matrix multiplication, second perspective
C1 = A[:,[0]]@B[[0],:]
C2 = A[:,[1]]@B[[1],:]
C3 = A[:,[2]]@B[[2],:]
fig, axs = plt.subplots(1, 7, figsize=(12, 3))
plt.sca(axs[0])
ax = sns.heatmap(C1,cmap='RdBu_r',vmax = all_max,vmin = all_min,
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,m+1), xticklabels=np.arange(1,n+1))
ax.set_aspect("equal")
plt.title('$C_1$')
plt.yticks(rotation=0)
plt.sca(axs[1])
plt.title('$+$')
plt.axis('off')
plt.sca(axs[2])
ax = sns.heatmap(C2,cmap='RdBu_r',vmax = all_max,vmin = all_min,
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,m+1), xticklabels=np.arange(1,n+1))
ax.set_aspect("equal")
plt.title('$C_2$')
plt.yticks(rotation=0)
plt.sca(axs[3])
plt.title('$+$')
plt.axis('off')
plt.sca(axs[4])
ax = sns.heatmap(C3,cmap='RdBu_r',vmax = all_max,vmin = all_min,
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,m+1), xticklabels=np.arange(1,n+1))
ax.set_aspect("equal")
plt.title('$C_3$')
plt.yticks(rotation=0)
plt.sca(axs[5])
plt.title('$=$')
plt.axis('off')
plt.sca(axs[6])
ax = sns.heatmap(C,cmap='RdBu_r',vmax = all_max,vmin = all_min,
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,m+1), xticklabels=np.arange(1,n+1))
ax.set_aspect("equal")
plt.title('$C$')
plt.yticks(rotation=0)
plt.show()
矩阵除法：计算逆矩阵
我们很容易想到在一个计算系统中，我们有单位元和逆元。对于矩阵的逆元运算，就是矩阵求逆。
逆运算的一些法则
广义逆矩阵
说明使用inv求逆的时候，此处能实现的是针对的可逆的矩阵，是方阵，且行列式不为零。
对格式进行检查
import numpy as np
A = np.array([[1, 2], [3, 4]])
det = np.linalg.det(A)
if det != 0:
A_inv = np.linalg.inv(A)
print("Inverse matrix:\n", A_inv)
else:
print("Matrix is singular and cannot be inverted.")
此外对于大规模的矩阵或者数值精度要求高的使用伪逆矩阵（广义逆矩阵）
补充：伪逆矩阵（pseudoinverse）就是广义逆矩阵的一个常用形式，通常指的是Moore-Penrose广义逆矩阵。它是用于解决那些不可逆或非方矩阵的逆问题的一种矩阵。
伪逆矩阵的主要特点是：
- 可以用于非方矩阵，即矩阵的行数和列数不相等的情况。
- 对于奇异矩阵（行列式为0的方阵），即便无法直接求逆，伪逆矩阵仍能给出一个广义的逆解。
在数值计算中，伪逆矩阵最常用于解决最小二乘法问题和处理欠定或超定系统。
在 numpy 中，np.linalg.pinv 就是计算Moore-Penrose伪逆的函数。其使用方法如下：
import numpy as np
A = np.array([[1, 2], [3, 4], [5, 6]]) # 非方矩阵
A_pinv = np.linalg.pinv(A)
print("Pseudoinverse matrix:\n", A_pinv)
Moore-Penrose广义逆矩阵满足以下四个条件：
- A A + A = A A A^+ A = A AA+A=A （矩阵乘积返回原矩阵）
- A + A A + = A + A^+ A A^+ = A^+ A+AA+=A+ （伪逆乘积返回伪逆矩阵）
- ( A A + ) T = A A + (A A^+)^T = A A^+ (AA+)T=AA+ （对称性条件）
- ( A + A ) T = A + A (A^+ A)^T = A^+ A (A+A)T=A+A （对称性条件）
因此，伪逆矩阵可以用来处理普通逆矩阵无法解决的情况，比如当矩阵不满足可逆性要求时。
代码inv求非奇异方阵的逆
#导入库
import numpy as np
#矩阵的逆运算
from numpy.linalg import inv
from matplotlib import pyplot as plt
import seaborn as sns
# Repeatability
np.random.seed(0)
# Generate matrix A
n = 4
A = np.random.uniform(-1.5,1.5,n*n).reshape(n, n)
all_max = 1.5
all_min = -1.5
# matrix inverse
A_inverse = inv(A)
#五个图像，三个矩阵，两个符号
fig, axs = plt.subplots(1, 5, figsize=(12, 3))
#这一行设置了当前的坐标轴为figure中的第一个子图（axs[0]）。在处理子图时这是必要的，因为每个子图都需要其自己的坐标轴
plt.sca(axs[0])
#使用seaborn库的heatmap函数创建热力图。A矩阵作为第一个参数传递，并传递了其他一些关键字参数：
#cmap='RdBu_r'：用于热力图的色图。'RdBu_r'是一个红蓝的色图，适用于表示具有正负值的数据。
#vmax = all_max：色图的最大值。在这里，它被设置为A矩阵中的最大值。
#vmin = all_min：色图的最小值。在这里，它被设置为A矩阵中的最小值。
#cbar_kws={"orientation": "horizontal"}：这一行设置了色图的水平和垂直方向。
#yticklabels=np.arange(1,n+1), xticklabels=np.arange(1,n+1),：这两行设置了y轴和x轴的刻度标签为从1到n+1的范围，其中n是A矩阵的大小。
#annot = True,fmt=".2f"：这两行启用了热力图单元格的注释，并设置了注释的格式为2位小数。
#ax.set_aspect("equal")：这一行设置了热力图的等高，即每个单元格都应该被一个正方形代表。
ax = sns.heatmap(A,cmap='RdBu_r',vmax = all_max,vmin = all_min,
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,n+1),
xticklabels=np.arange(1,n+1),
annot = True,fmt=".2f")
ax.set_aspect("equal")
plt.title('$A$')
plt.yticks(rotation=0)
#设置符号
plt.sca(axs[1])
plt.title('$@$')
plt.axis('off')
plt.sca(axs[2])
ax = sns.heatmap(A_inverse,cmap='RdBu_r',vmax = all_max,vmin = all_min,
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,n+1), xticklabels=np.arange(1,n+1),
annot = True,fmt=".2f")
ax.set_aspect("equal")
plt.title('$A^{-1}$')
plt.yticks(rotation=0)
plt.sca(axs[3])
plt.title('$=$')
plt.axis('off')
plt.sca(axs[4])
ax = sns.heatmap(A@A_inverse,cmap='RdBu_r',vmax = all_max,vmin = all_min,
cbar_kws={"orientation": "horizontal"},
yticklabels=np.arange(1,n+1), xticklabels=np.arange(1,n+1),
annot = True,fmt=".2f")
ax.set_aspect("equal")
plt.title('$I$')
plt.yticks(rotation=0)
plt.show()
