python中是否有一种有效的方法将矩阵中的每一列与同一矩阵中的所有列相乘？

Question

我有一个大矩阵（例如 100.000 x 100.000）。好在它只包含零和一，而且大部分都是零（它已经保存为布尔矩阵以节省一些 RAM）。现在我需要将矩阵的每一列与所有其他列相乘。原因是我需要检查是否至少有一行两列都具有非零元素（因此将结果向量相乘并求和以检查它是否为零）。例如假设我们有一个矩阵

1.column	2.column	3.column
1	0	0
1	1	0
0	0	1

然后我需要比较所有列并检查是否至少有一行两列都为一。因此，比较第一列和第二列将返回 True，因为它们都是第二行中的一个。但是，比较第一列和第三列以及第二列第三列将导致 False，因为没有行的行都为一个。 显然，这可以使用 for 循环并遍历所有列来完成。但是速度不是很令人满意。我已经尝试过这样的 numba：

@njit(parallel=True)
def create_dist_arr(arr: np.array):
    n = arr.shape[1]
    dist_arr = np.zeros(shape=(n, n)) #, dtype=bool)
    for i in prange(arr.shape[1]):
        for j in prange(i, arr.shape[1]):
            dist_greater_zero = calc_dist_graeter_than_zeros(arr[:, i], arr[:, j])
            dist_arr[i][j] = dist_greater_zero
            dist_arr[i][j] = dist_greater_zero
    return skill_dist_arr

@njit
def calc_dist_graeter_than_zeros(ith_col, jth_col):
    return np.sum(np.multiply(ith_col, jth_col)) != 0

zero_arr = np.zeros(shape=(2000, 6000), dtype=bool)
bool_dist_matrix = create_dist_arr(zero_arr)

但是尽管有 120gb 的 RAM 和 32 个内核，但在 10.000 x 10.000 矩阵左右时会变得非常慢。像这样尝试 scipy.spatial.distance.pdist 时更糟糕的是：

from scipy.spatial.distance import pdist
zero_arr = np.zeros(shape=(500, 500), dtype=bool)
bool_dist_matrix = pdist(zero_arr, lambda u, v: np.sum(np.multiply(u, v)) != 0)

是否有一些使用稀疏矩阵或其他不会永远使用的好而快速的解决方法？

提前谢谢你:)

Answer 1

想知道这个解决方案的内存效率如何，也不知道它是否比其他解决方案更快，但它是矢量化的。

您的想法是将列相乘并添加结果。这让我想起了矩阵乘法，除了它针对自己的列。所以..

假设您的矩阵是M。如果您将横向 M.T 和矩阵相乘，M.T @ M 将与将每一列与其他列相乘并求和相同。

import numpy as np

# A random matrix
M = np.array([[1,0,0,0,0,1],
              [1,1,0,0,0,1],
              [0,1,1,1,0,0]])

bool_dist_matrix = (M.T @ M).astype('bool')
np.fill_diagonal(bool_dist_matrix, 0)

"""
[[0 1 0 0 0 1]
 [1 0 1 1 0 1]
 [0 1 0 1 0 0]
 [0 1 1 0 0 0]
 [0 0 0 0 0 0]
 [1 1 0 0 0 0]]
"""

Answer 2

如果你的矩阵主要由零组成，那么使用索引会更有效：

import numpy as np

# A random matrix
M = np.array([[1,0,0,0,0,1],
              [1,1,0,0,0,1],
              [0,1,1,1,0,0]])

# Get the index where M == 1
ind = np.where(M)
# Get the unique value and return the count.
uni = np.unique(ind[0], return_counts=True)
# Keep only the column that have at least two 1 in the same row.
col = uni[0][uni[1]>1]
# Create a boolean index.
bol = np.isin(ind[0],col)

# And here we get the "compressed" result:
res = ind[1][bol] #Col number
grp = ind[0][bol] #Group
# res = [0, 5, 0, 1, 5, 1, 2, 3]
# grp = [1, 1, 2, 2, 2, 3, 3, 3]
# So each pair of each group will output a True statement:
# group 1: 0-5
# group 2: 0-1, 0-5, 1-5
# group 3: 1-2, 1-3, 2-3
# For an explicit result you could use itertools. But all the information is 
# contained in the res variable.

注意到，如果对于多行，两列具有共同的 1 值，则此方法会产生一些重复。但是很容易摆脱那些重复的。但由于您使用的是 100000x100000 矩阵，并且并非所有列都至少有一个共同的 1 值，因此矩阵中 1 的百分比很可能非常低，因此这种方法应该会提供一些不错的结果。

Answer 3

我不确定你的意思是你需要做什么，但如果我理解正确，我的代码应该会有所帮助。这是我尝试过的。它似乎跑得快了一点。虽然如果矩阵不像你提到的那样稀疏并且我得到一个对称矩阵而不是上三角矩阵，则情况并非如此：

import numpy as np
from numba import njit, prange
@njit(parallel=True)
def create_dist_arr(arr: np.array):
    n = arr.shape[1]
    dist_arr = np.zeros(shape=(n, n)) #, dtype=bool)
    for i in prange(arr.shape[1]):
        for j in prange(i, arr.shape[1]):
            dist_greater_zero = calc_dist_graeter_than_zeros(arr[:, i], arr[:, j])
            dist_arr[i][j] = dist_greater_zero
            dist_arr[i][j] = dist_greater_zero
    return dist_arr

@njit
def calc_dist_graeter_than_zeros(ith_col, jth_col):
    return np.sum(np.multiply(ith_col, jth_col)) != 0

def create_dist_arr_sparse(arr: np.array):
    n = arr.shape[1]
    dist_arr = np.zeros(shape=(n, n)) #, dtype=bool)
    rows, cols = np.array(np.where(arr))
    same_rows = rows.reshape(1, -1) == rows.reshape(-1, 1)
    idx0, idx1 = np.meshgrid(cols, cols)
    idx0 = idx0.flatten()[same_rows.flatten()]
    idx1 = idx1.flatten()[same_rows.flatten()]
    dist_arr[idx0, idx1] = 1
    return dist_arr

np.random.seed(1)
k = 1000
zero_arr = np.zeros(shape=(k, k), dtype=bool)
rows, cols = np.random.randint(0, k, (2, k))
zero_arr[rows, cols] = 1
%timeit bool_dist_matrix = create_dist_arr(zero_arr)
%timeit bool_dist_matrix = create_dist_arr_sparse(zero_arr)

输出：

1 loop, best of 5: 1.59 s per loop
100 loops, best of 5: 9.8 ms per loop