Python，Scipy：使用大型邻接矩阵构建三元组

Question

我正在使用邻接矩阵来表示一个朋友网络，可以被视为

Mary     0        1      1      1

Joe      1        0      1      1

Bob      1        1      0      1

Susan    1        1      1      0 

         Mary     Joe    Bob    Susan

我想使用这个矩阵编译出所有可能的友谊三角形列表，条件是用户1和用户2是朋友，用户2和用户3也是朋友。对于我的列表，不需要要求用户1和用户3是朋友。

(joe, mary, bob)
(joe, mary, susan)
(bob, mary, susan)
(bob, joe, susan)

我有一段代码适用于小三角形，但需要对非常大的稀疏矩阵进行扩展。

from numpy import *
from scipy import *

def buildTriangles(G):
    # G is a sparse adjacency matrix
    start = time.time()
    ctr = 0
    G = G + G.T          # I do this to make sure it is symmetric
    triples = []
    for i in arange(G.shape[0] - 1):  # for each row but the last one
        J,J = G[i,:].nonzero()        # J: primary friends of user i
                                      # I do J,J because I do not care about the row values
        J = J[ J < i ]                # only computer the lower triangle to avoid repetition
        for j in J:
            K, buff = G[:,j].nonzero() # K: secondary friends of user i
            K = K[ K > i ]             # only compute below i to avoid repetition
            for k in K:
                ctr = ctr + 1
                triples.append( (i,j,k) )
    print("total number of triples: %d" % ctr)
    print("run time is %.2f" % (time.time() - start())
    return triples

我能够在约21分钟内对1032570 x 1032570的csr_matrix运行代码，其中包含88910个存储元素。总共生成了2178893个三元组。

我需要能够处理类似于1968654 x 1968654的稀疏矩阵，其中包含9428596个存储元素。

我对Python非常陌生（不到一个月的经验），而且不是线性代数方面的专家，这就是为什么我的代码没有利用矩阵操作的原因。有人能提出改进建议或让我知道我的目标是否可行吗？

- will

---

我认为在Python中在一个语句中两次赋相同的值（J，J =）没有任何保证的意义。我发现这非常令人困惑，你的评论也是如此，所以你可能想要摆脱它。 - Fred Foo

@larsmans 抱歉。nonzero()函数返回矩阵的索引作为2维数组。或者我可以这样做：row，col = G[i,:] .nonzero()，然后J = col。我使用J，J=方法是因为我担心内存使用情况，并希望消耗掉行数组，因为它不需要。 - will

不用道歉，我并不是要刻意苛刻。只是这不是 Python 惯用法，并且我认为 Guido 有权更改该构造在 Python 版本之间的含义，因此不能依赖它的工作方式。如果确实很重要，最好使用 del 删除变量，尽管在这种情况下，J = G[i, :].nonzero()[1] 也可以工作。 - Fred Foo

谢谢您的建议。它确实使代码变得更简洁了。您正在处理维基百科文章的工作正是我想要做的。我将进一步研究线性代数方法来解决这个问题。 - will

Answer 1

我认为你只能在行或列中找到三角形，例如：

Susan    1        1      1      0 
        Mary     Joe    Bob    Susan

这意味着玛丽、乔和鲍勃都是苏珊的朋友，因此，使用组合从[Mary，Joe，Bob]中选择两个人，并将其与Susan结合起来，将得到一个三角形。 itertools.combinations() 可以快速实现这一点。

以下是代码：

import itertools
import numpy as np

G = np.array(   # clear half of the matrix first
    [[0,0,0,0],
     [1,0,0,0],
     [1,1,0,0],
     [1,1,1,0]])
triples = []     
for i in xrange(G.shape[0]):
    row = G[i,:]
    J = np.nonzero(row)[0].tolist() # combinations() with list is faster than NumPy array.
    for t1,t2 in itertools.combinations(J, 2):
        triples.append((i,t1,t2))
print triples

- HYRY

---

谢谢你的回答。我甚至没有考虑过那种方法，但它非常有道理。你基本上将问题简化为找到两个的排列。所有的三元组都是唯一的吗？ - will

@will：请澄清一下，您指的是(Mary, Susan, Joe)和(Joe, Susan, Mary)被视为不同还是相同？ - Iterator

@Iterator 我的意思是将它们视为相同。我相信这种方法在这方面确实有效。经过进一步的研究，我现在意识到每一行都保证不会出现在早期的排列中。 - will

+1 给用户772649。这太棒了。我想在我工作的其他语言中找到这个函数。我总是不得不自己编写它。 - Iterator

Answer 2

以下是优化建议：

以下是一些优化建议：

K = K[ K > i ]             # only compute below i to avoid repetition
for k in K:
    ctr = ctr + 1
    triples.append( (i,j,k) )

不要在循环中递增，这会使速度变得非常慢。只需要使用ctr += K.shape[0]即可。然后，通过将append替换为...来完全消除最深层的循环。

triples += ((i, j, k) for k in K[K > i])

现在，如果你想要在这个任务上获得真正的性能，你需要涉及一些线性代数。"我想编制一个可能的友谊三角形列表"意味着你想要平方邻接矩阵，可以通过简单的**2来完成。

然后请注意，1,968,654²意味着一个非常大的矩阵，即使它非常稀疏，它的平方也会少得多，并且需要大量的内存。(我曾经处理过一个类似的问题，在其中考虑了维基百科文章之间的距离二的链接，在C++中，在超级计算机集群节点上解决它需要20分钟。这不是一个简单的问题。尽管维基百科邻接矩阵密度高几个数量级。)

- Fred Foo

---

当您提到“真实性能”时，您能否详细说明如何将两个矩阵相乘并获得2步配对的列表（而不是计数）？ - Iterator

@迭代器：将一个方阵与其自身相乘，可以得到一个新的同秩矩阵，其中所有连接在步长2处的i和j都具有值>0。矩阵乘法是SciPy中高度优化的操作（我认为是用C实现的，甚至可能是Fortran）。然后，您可以通过在矩阵中进行更少的搜索来提取列表。 - Fred Foo

是的，你得到了第二步的计数，也就是我说的：你可以得到(i,*,k)对的计数。中间j节点的身份丢失了。我理解（并声明）你所说的一切，但你没有展示出完整三元组的命名加速。我认为你没有完全思考清楚这个问题。 - Iterator