Python求矩阵的内积、外积、克罗内克直积、Khatri-Rao积

矩阵乘法

线性代数研究的核心对象是矩阵，所谓矩阵就是由$m$行$n$列的数组成的一个举行的数阵，从编程的角度理解，就是二维数组。

作为一种特殊的代数结构，需要为矩阵引入特殊的计算方法，其中矩阵加法、减法相对来说比较直观，而乘法比较复杂，矩阵乘法的定义为

(AB)ij=∑aikbkj(AB)_{ij}=\\sum a_{ik}b_{kj} (AB)ij​=∑aik​bkj​

在python中可直接通过$A@B$来实现，而常用的乘号*表示元素之间一对一的乘法。

此外，numpy中提供了np.matmul进行矩阵乘法。

内积和外积

对于向量$a,b$，其内积表示为

a⋅b=∑iaibia\\cdot b=\\sum_i a_ib_i a⋅b=i∑​ai​bi​

外积则得到一个矩阵

(AB)ij=aibj(AB)_{ij}=a_ib_j (AB)ij​=ai​bj​

而矩阵乘法则可表示为向量外积之和。

在numpy中，提供了内积inner和外积outer，令a=[1234],b=[0321]a=\\begin{bmatrix}1&2\\\\3&4\\end{bmatrix}, b=\\begin{bmatrix}0&3\\\\2&1\\end{bmatrix}a=[13​24​],b=[02​31​]，则计算结果如下

import numpy as np
a = np.array([1,2,3,1]).reshape(2,2)
b = np.array([0,3,2,1]).reshape(2,2)
print(np.inner(a,b))
'''
[[6 4][3 7]]
'''
print(np.outer(a,b))
'''
[[0 3 2 1][0 6 4 2][0 9 6 3][0 3 2 1]]
'''

直积

numpy和scipy.linalg中都实现了矩阵直积，又称克罗内克积

A⊗B=[a11Ba12B⋯a1nBa21Ba22B⋯a2nB⋯⋯⋯⋯am1Bam2B⋯amnB]A\\otimes B=\\begin{bmatrix} a_{11}B&a_{12}B&\\cdots&a_{1n}B\\\\ a_{21}B&a_{22}B&\\cdots&a_{2n}B\\\\ \\cdots&\\cdots&\\cdots&\\cdots\\\\ a_{m1}B&a_{m2}B&\\cdots&a_{mn}B\\\\ \\end{bmatrix} A⊗B=​a11​Ba21​B⋯am1​B​a12​Ba22​B⋯am2​B​⋯⋯⋯⋯​a1n​Ba2n​B⋯amn​B​​

import scipy.linalg as sl
c = sl.kron(a,b)
print(c)
'''
[[0 3 0 6][2 1 4 2][0 9 0 3][6 3 2 1]]
'''
print(np.kron(a,b))
'''
[[0 3 0 6][2 1 4 2][0 9 0 3][6 3 2 1]]
'''

下面来看一下二者的速度

>>> import timeit
>>> y = np.random.rand(30,30)
>>> x = np.random.rand(30,30)
>>> timeit.timeit(lambda :np.kron(x,y), number=10)
0.023692000002483837
>>> timeit.timeit(lambda :sl.kron(x,y), number=10)
0.05632819999300409

看来这次是scipy输了。

Khatri-Rao积

Khatri-Rao积是两个具有相同列数的矩阵$A,B$对应列向量的克罗内克积排列而成的，可表示为

A⊙B=[a1⊗b1a2⊗b2⋯an⊗bn]A\\odot B = \\begin{bmatrix} a_1\\otimes b_1&a_2\\otimes b_2&\\cdots&a_n\\otimes b_n \\end{bmatrix} A⊙B=[a1​⊗b1​​a2​⊗b2​​⋯​an​⊗bn​​]

Khatri-Rao积满足结合律，在scipy.linalg中，用khatri_rao可求矩阵的Khatri-Rao积。

d = sl.khatri_rao(a,b)
print(d)
'''
[[0 6][2 2][0 3][6 1]]
'''