稀疏矩阵之 toarray 方法和 todense 方法

不可言诉的深渊

发布于 2022-11-11 15:02:23

3.3K0

发布于 2022-11-11 15:02:23

在 SciPy 稀疏矩阵中，有着 2 个经常被混为一谈的方法：toarray() 方法以及 todense() 方法。事实上，我在才开始接触 SciPy 稀疏矩阵的时候也曾经把这 2 个方法之间画上等号。但是，两者之间还是存在着很大的不同，具体有哪些不同之处我们就首先从返回值类型开始说明。

返回值类型

在说明返回值类型之前，我们首先需要知道的是不管是 toarray() 方法还是 todense() 方法，它们都是 7 种 SciPy 稀疏矩阵中的任意一种稀疏矩阵类的实例的方法！接下来我就以 COO 格式的 SciPy 稀疏矩阵作为示例说明一下 toarray() 方法以及 todense() 方法的返回值，代码如下：

>>> import numpy as np
>>> from scipy.sparse import coo_matrix
>>> random_state = np.random.RandomState(0)
>>> row = random_state.randint(8, size=8)
>>> col = random_state.randint(8, size=8)
>>> data = random_state.randint(-4, 4, 8)
>>> coo = coo_matrix((data, (row, col)), shape=(8, 8))
>>> a = coo.toarray()
>>> b = coo.todense()
>>> a
array([[ 0,  0, -3,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  2,  0,  3,  3],
       [ 0, -4,  0,  0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0, -2,  0,  0],
       [ 0,  0,  0,  0,  0,  0,  0,  0],
       [ 2,  0,  0,  0,  0,  0,  0,  0]])
>>> b
matrix([[ 0,  0, -3,  0,  0,  0,  0,  0],
        [ 0,  0,  0,  0,  0,  0,  0,  0],
        [ 0,  0,  0,  0,  0,  0,  0,  0],
        [ 0,  0,  0,  0,  2,  0,  3,  3],
        [ 0, -4,  0,  0,  0,  0,  0,  0],
        [ 0,  0,  0,  0,  0, -2,  0,  0],
        [ 0,  0,  0,  0,  0,  0,  0,  0],
        [ 2,  0,  0,  0,  0,  0,  0,  0]])
>>> type(a)
<class 'numpy.ndarray'>
>>> type(b)
<class 'numpy.matrix'>

需要注意的是 Python 和 C 语言不一样，定义函数的时候完全不需要指定返回值的类型，调用函数的时候接收返回值的变量也同样是完全不需要指定其对应的类型。这应该大概可能也许就是让 SciPy 稀疏矩阵的初学者把二者混为一谈的主要原因吧。显然，我们可以发现 toarray() 方法的返回值类型是 numpy.ndarray 而 todense() 方法的返回值类型是 numpy.matrix，这两个类型必然存在某种程度上的不同之处。至于为什么我们可以反过来想一下，如果这两个类型完全一样的话，那么 NumPy 内部要去实现这两个类型就会显得非常的冗余。因此，为了避开所谓的冗余的代码，NumPy 内部对两者的内部实现存在着某种程度上的不同。

numpy.ndarray：numpy 就是模块名，我们完全没有必要去管它。ndarray 要拆成 3 个部分，n 表示 Natural Number（自然数，其值我们就简单的记作变量 n），d 表示Dimension（维度），至于 array 嘛，就不用我多说了吧，直接翻译成中文即可，它就是数组。因此，numpy.ndarray 表示 NumPy 模块中的 n 维数组类。

numpy.matrix：numpy 和上面一样，也是模块名，我们依旧完全没有必要去管它。至于 matrix，它没有上面那么复杂，和 array 一样，直接翻译成中文就行了，它就是矩阵。因此，numpy.matrix 表示 NumPy 模块中的矩阵类。

因此，当 ndarray 的 n≠2 的时候，ndarray 类的某个实例就绝对不可能是一个矩阵，至少无法看作是一个矩阵。接下来我需要重点说明的是当 ndarray 的 n=2 的时候，它和 matrix 又存在着怎么样的区别。

二维数组和矩阵

当 ndarray 的 n=2 的时候，我们可以把这个类的某个实例当成是一个矩阵，因此，我们猜测可以对其进行矩阵运算。我们都知道矩阵的运算无非就是加法、减法、数乘、转置、乘法、求逆、求幂、哈达玛乘积和克罗内克乘积。其中，加法、减法、乘法、哈达玛乘积和克罗内克乘积是二元运算，两个操作变量都是矩阵；数乘运算也是二元运算，只不过它的两个操作变量是一个数和一个矩阵；转置、求逆和求幂都是一元运算，操作变量只有一个矩阵。在这些运算中，我们需要注意的是加法、减法和哈达玛乘积必须确保两个矩阵形状相同；乘法运算必须确保第一个矩阵的列数和第二个矩阵的行数必须完全相等；求逆运算必须确保矩阵是一个可逆方阵；求幂运算，求的是方阵的幂，我们假设求的是 n 次幂，n 是一个整数（可正可负）。当 n<0 的时候，返回把 |n| 个对应方阵的逆矩阵用矩阵乘法连接在一起进行运算得到的结果；当 n=0 的时候，返回和对应方阵形状相同的单位矩阵（主对角线全 1，其余位置一律为 0）；当 n>0 的时候，返回把 n 个对应方阵用矩阵乘法连接在一起进行运算得到的结果。因此，当 n≥0 的时候，仅仅要求该矩阵是一个方阵就行了；而当 n<0 的时候，它不仅要求该矩阵是一个方阵，而且还加上了可逆的限制条件。至于数乘、转置和克罗内克乘积，它们 3 个运算对矩阵的要求没有任何限制。

二维数组

我们首先看一下把上述的几种运算丢给二维数组看看会得出什么样的结果，代码及结果如下所示：

>>> import numpy as np
>>> random_state = np.random.RandomState(0)
>>> a = random_state.rand(2, 2)
>>> a
array([[0.5488135 , 0.71518937],
       [0.60276338, 0.54488318]])
>>> a+a
array([[1.09762701, 1.43037873],
       [1.20552675, 1.08976637]])
>>> a-a
array([[0., 0.],
       [0., 0.]])
>>> 2*a
array([[1.09762701, 1.43037873],
       [1.20552675, 1.08976637]])
>>> a.T
array([[0.5488135 , 0.60276338],
       [0.71518937, 0.54488318]])
>>> a*a
array([[0.30119626, 0.51149583],
       [0.36332369, 0.29689768]])
>>> a**-1
array([[1.8221126 , 1.39823108],
       [1.65902581, 1.83525576]])
>>> a**2
array([[0.30119626, 0.51149583],
       [0.36332369, 0.29689768]])
>>> a@a
array([[0.73228622, 0.78220024],
       [0.65924031, 0.72798764]])
>>> np.linalg.inv(a)
array([[-4.12631777,  5.41602068],
       [ 4.5646357 , -4.15608149]])
>>> np.kron(a, a)
array([[0.30119626, 0.39250558, 0.39250558, 0.51149583],
       [0.33080468, 0.29903925, 0.43108996, 0.38969466],
       [0.33080468, 0.43108996, 0.29903925, 0.38969466],
       [0.36332369, 0.32843563, 0.32843563, 0.29689768]])

显然我们可以发现二维数组的加法相当于矩阵的加法，二维数组的减法相当于矩阵的减法，一个数乘上一个二维数组相当于一个数乘上一个矩阵，二维数组的转置相当于矩阵的转置。但是，第一，二维数组的乘法和矩阵的乘法并不能划等号，二维数组的乘法是把两个相同形状的二维数组的对应位置的元素相乘得到一个新数组，和矩阵的乘法并不能画上等号，如果把二维数组看作是矩阵，这就相当于两个矩阵做哈达玛乘积；第二，二维数组的 -1 次方和矩阵的逆也不能画上等号，二维数组的 -1 次方是在对二维数组中的每个元素计算 -1 次方并得到一个新的二维数组；第三，二维数组的 n 次幂也同样不等于矩阵的 n 次幂，二维数组的 n 次幂是在对二维数组中的每个元素计算 n 次幂并得到一个新的二维数组。我们可以针对它来实现对应的矩阵乘法、矩阵的逆以及矩阵的克罗内克乘积，矩阵乘法很简单，把 * 运算符改成 @ 运算符就行了；矩阵的逆就需要调用 np.linalg.inv 函数，参数就是需要求逆的矩阵（二维数组）；矩阵的克罗内克乘积需要调用 np.kron 函数，第一个参数是左操作变量，第二个参数是右操作变量。至于针对二维数组来计算其对应矩阵的幂，我没有找到对应的函数，但是我们可以借助上面的 @ 运算来手工实现，代码如下：

>>> def matrix_power(a, n):
...     assert a.ndim == 2 and a.shape[0] == a.shape[1] and isinstance(n, int)
...     if n > 0:
...         b = a.copy()
...         for _ in range(n-1):
...             b = b@a
...     elif n < 0:
...         b = inv_a = np.linalg.inv(a)
...         for _ in range(-n-1):
...             b = b@inv_a
...     else:
...         b = np.eye(a.shape[0])
...     return b
...
>>> matrix_power(a, 3)
array([[0.87337022, 0.94993107],
       [0.80060427, 0.86814988]])
>>> a@a@a
array([[0.87337022, 0.94993107],
       [0.80060427, 0.86814988]])
>>> inv_a = np.linalg.inv(a)
>>> matrix_power(a, -3)
array([[-377.02704688,  412.5436352 ],
       [ 347.69280141, -379.29417914]])
>>> inv_a@inv_a@inv_a
array([[-377.02704688,  412.5436352 ],
       [ 347.69280141, -379.29417914]])
>>> matrix_power(a, 0)
array([[1., 0.],
       [0., 1.]])

‍

首先说一下函数 matrix_power 的具体内容，该方法用来求二维数组 a 所对应矩阵的 n 次幂。我来详细说明一下这个函数，首先是该函数的两个参数，参数 a 是一个 numpy.ndarray 类的实例，参数 n 是一个整数。然后看到函数体的第 1 行的 assert 断言，这里通过逻辑与 and 运算去连接 3 个条件，第 1 个条件用来判断 a 是不是二维数组（换句话说，它用来判断 a 是不是可以被看作是一个矩阵），第 2 个条件用来判断 a 所对应的矩阵是不是方阵，第 3 个条件用来判断参数 n 是不是整数。当这 3 个条件都为真的时候才能进行矩阵的求幂运算。求矩阵的 n 次幂需要分成 3 种情况：n 为正整数，n 为负整数，n=0。这也就对应着函数体内的 3 个互斥条件分支。当 n≠0 的时候，它要把 |n| 个 a（或者 a 的逆）通过矩阵乘法连在一起，也就是矩阵连乘 |n|-1 次，所以不管是 n>0 还是 n<0 所对应的条件分支都是循环 |n|-1 次就够了。需要注意的是，第一，当 n<0 的时候，我之所以没有去判断矩阵 a 是否可逆是因为在 np.linalg.inv 函数内部已经实现了类似的判断；第二，第一个条件分支里面之所以是 b = a.copy() 而不是 b = a 是因为我不希望返回值和参数 a 捆绑在一起（共享同样的内存空间，当 n=1 的时候就会出现这种情况）。

通过观察针对该函数的简单测试，我们可以发现它可以在不修改参数类型的情况下实现二维数组所对应矩阵的 n 次幂，二维数组自始至终都是二维数组，没有转换为矩阵（numpy.matrix 类的实例）。

关于二维数组及其对应矩阵的各种运算就介绍到这里，接下来我们来看矩阵类（numpy.matrix）的实例所对应的各种运算。

矩阵

在讲矩阵运算之前，我们首先需要看一下通过一个二维数组来构造一个矩阵的方法，这样的方法有很多，我比较推荐去使用 numpy.mat 函数，这个函数接受一个参数，该参数就是二维数组。对矩阵进行之前提到的运算代码和结果如下所示：

>>> import numpy as np
>>> random_state = np.random.RandomState(0)
>>> a = np.mat(random_state.rand(2, 2))
>>> a
matrix([[0.5488135 , 0.71518937],
        [0.60276338, 0.54488318]])
>>> a+a
matrix([[1.09762701, 1.43037873],
        [1.20552675, 1.08976637]])
>>> a-a
matrix([[0., 0.],
        [0., 0.]])
>>> 2*a
matrix([[1.09762701, 1.43037873],
        [1.20552675, 1.08976637]])
>>> a.T
matrix([[0.5488135 , 0.60276338],
        [0.71518937, 0.54488318]])
>>> a*a
matrix([[0.73228622, 0.78220024],
        [0.65924031, 0.72798764]])
>>> a**-1
matrix([[-4.12631777,  5.41602068],
        [ 4.5646357 , -4.15608149]])
>>> a**2
matrix([[0.73228622, 0.78220024],
        [0.65924031, 0.72798764]])
>>> a@a
matrix([[0.73228622, 0.78220024],
        [0.65924031, 0.72798764]])
>>> np.multiply(a, a)
matrix([[0.30119626, 0.51149583],
        [0.36332369, 0.29689768]])
>>> np.power(a, 2)
matrix([[0.30119626, 0.51149583],
        [0.36332369, 0.29689768]])
>>> np.kron(a, a)
matrix([[0.30119626, 0.39250558, 0.39250558, 0.51149583],
        [0.33080468, 0.29903925, 0.43108996, 0.38969466],
        [0.33080468, 0.43108996, 0.29903925, 0.38969466],
        [0.36332369, 0.32843563, 0.32843563, 0.29689768]])

对于 numpy.matrix 类的实例来说，加法、减法、数乘以及转置全部都和二维数组的结果几乎完全一致。需要注意的是，第一，若 * 运算符的两个操作变量是两个矩阵，则它再也不是用来表示两个矩阵的哈达玛乘积，而是用来表示原生的矩阵乘法，和 @ 运算符完全等价；第二，若 ** 运算符的左操作变量是矩阵右操作变量是一个整数 n，则它再也不是用来表示矩阵中的每个元素求 n 次幂得到新矩阵，而是用来表示矩阵的原生的 n 次幂，当 n=-1 时求的就是矩阵的逆。如果我想要给矩阵实现二维数组的乘法（对应元素相乘，哈达玛乘积），那么可以调用 np.multiply 函数，两个参数类型都是矩阵；如果我想要给矩阵实现二维数组的 n 次幂，可以调用 np.power 函数，第一个参数是矩阵，第二个参数是整数 n；克罗内克乘积和二维数组一样，依旧是调用 np.kron 函数，两个参数类型都是矩阵。

混合运算

最后我们需要看一下如果两个操作变量其中一个是二维数组（numpy.ndarray 类的实例），而另一个是矩阵（numpy.matrix 类的实例），让它们进行之前提到的两个操作变量都是矩阵的二元运算会出现什么样的结果。代码和结果如下所示：

>>> import numpy as np
>>> random_state = np.random.RandomState(0)
>>> a = random_state.rand(2, 2)
>>> b = np.mat(random_state.rand(2, 2))
>>> a
array([[0.5488135 , 0.71518937],
       [0.60276338, 0.54488318]])
>>> b
matrix([[0.4236548 , 0.64589411],
        [0.43758721, 0.891773  ]])
>>> a+b
matrix([[0.9724683 , 1.36108348],
        [1.04035059, 1.43665618]])
>>> b+a
matrix([[0.9724683 , 1.36108348],
        [1.04035059, 1.43665618]])
>>> a-b
matrix([[ 0.1251587 ,  0.06929525],
        [ 0.16517616, -0.34688982]])
>>> b-a
matrix([[-0.1251587 , -0.06929525],
        [-0.16517616,  0.34688982]])
>>> a*b
matrix([[0.5454652 , 0.99226198],
        [0.49379751, 0.87523343]])
>>> b*a
matrix([[0.62182879, 0.65493025],
        [0.77768188, 0.79886983]])
>>> a@b
matrix([[0.5454652 , 0.99226198],
        [0.49379751, 0.87523343]])
>>> b@a
matrix([[0.62182879, 0.65493025],
        [0.77768188, 0.79886983]])
>>> np.multiply(a, b)
matrix([[0.23250747, 0.4619366 ],
        [0.26376154, 0.48591211]])
>>> np.multiply(b, a)
matrix([[0.23250747, 0.4619366 ],
        [0.26376154, 0.48591211]])
>>> np.kron(a, b)
matrix([[0.23250747, 0.35447541, 0.30299341, 0.4619366 ],
        [0.24015377, 0.48941707, 0.31295772, 0.63778657],
        [0.2553636 , 0.38932132, 0.23084238, 0.35193684],
        [0.26376154, 0.5375281 , 0.23843391, 0.48591211]])
>>> np.kron(b, a)
matrix([[0.23250747, 0.30299341, 0.35447541, 0.4619366 ],
        [0.2553636 , 0.23084238, 0.38932132, 0.35193684],
        [0.24015377, 0.31295772, 0.48941707, 0.63778657],
        [0.26376154, 0.23843391, 0.5375281 , 0.48591211]])

显然，第一，当存在一个矩阵一个二维数组的时候，在运算之前会先把二维数组给隐式地转换为矩阵再进行运算；第二，所有的运算符都按照矩阵的来。

结论

在这里，我首先通过稀疏矩阵的 toarray() 方法以及 todense() 方法的返回值看似一样但实际上却是两个完全不同的类的实例，然后通过对矩阵的运算给出它们两者的区别。最后给出一些注意事项：