比较两个不同大小的矩阵以制作一个大矩阵 - 速度改进？

Posted 2023-03-11

【中文标题】比较两个不同大小的矩阵以制作一个大矩阵 - 速度改进？

【英文标题】：Compare two different size matrices to make one large matrix - Speed Improvements?

【发布时间】：2016-10-03 18:28:44

【问题描述】：

我有两个矩阵，我需要用它们来创建一个更大的矩阵。每个矩阵只是一个以制表符分隔的文本文件，可供读取。每个矩阵有 48 个列，每个矩阵具有相同的标识符，具有不同的行数。第一个矩阵是 108887x48，第二个是 55482x48。对于每个矩阵，每个位置的条目可以是 0 或 1，因此是二进制的。最终输出应该有第一个矩阵的行 ID 作为行，第二个矩阵的行 ID 作为列，所以最终的矩阵是 55482x10887。

这里需要发生的是，对于第一个矩阵中每一行中的每个 pos，对于第二个矩阵中的每一行，如果每个矩阵的 pos(col) 为 1，那么最终的矩阵计数将上升 1 . 最终矩阵中任意 pos 的最大值为 48，预计会剩下 0。

例子：

mat1
     A B C D
1id1 0 1 0 1
1id2 1 1 0 0
1id3 1 1 1 1
1id4 0 0 1 0

mat2
     A B C D
2id1 1 1 0 0
2id2 0 1 1 0 
2id3 1 1 1 1 
2id4 1 0 1 0

final
     2id1 2id2 2id3 2id4
1id1   1    1    2    0
1id2   2    1    2    1
1id3   2    2    4    2
1id4   0    1    1    1

我有代码可以做到这一点，但是速度非常慢，这是我主要寻求帮助的地方。我试图尽可能地加速算法。它已经运行了 24 小时，仅完成了大约 25%。我之前已经让它运行过了，最终的输出文件是 20GB。我对数据库没有经验，可以在这里实现它，如果有人可以在给定下面代码的 sn-p 的情况下帮助我如何做到这一点。

#!/usr/bin/env python

import sys

mat1in = sys.argv[1]
mat2in = sys.argv[2]

print '\n######################################################################################'
print 'Generating matrix by counts from smaller matrices.'
print '########################################################################################\n'

with open(mat1in, 'r') as f:
        cols = [''] + next(f).strip().split('\t')               # First line of matrix is composed of 48 cols
        mat1 = [line.strip().split('\t') for line in f]         # Each line in matrix = 'ID': 0 or 1 per col id

with open(mat2in, 'r') as f:
        next(f)                                                 # Skip first row, col IDs are taken from mat1
        mat2 = [line.strip().split('\t') for line in f]         # Each line in matrix = 'ID': 0 or 1 per col id

out = open('final_matrix.txt', 'w')                             # Output file

#matrix = []
header = []                                                     # Final matrix header
header.append('')                                               # Add blank as first char in large matrix header
for i in mat2:
        header.append(i[0])                                     # Composed of all mat2 row ids
#matrix.append(header)

print >> out, '\t'.join(header)                                 # First print header to output file

print '\nTotal mat1 rows: ' + str(len(mat1))                    # Get total mat1 rows
print 'Total mat2 rows: ' + str(len(mat2)), '\n'                # Get total mat2 rows
print 'Progress: '                                              # Progress updated as each mat1 id is read

length = len(header)                                            # Length of header, i.e. total number of mat2 ids
totmat1 = len(mat1)                                             # Length of rows (-header), i.e. total number of mat1 ids

total = 0                                                       # Running total - for progress meter
for h in mat1:                                                  # Loop through all mat1 ids - each row in the HC matrix
        row = []                                                # Empty list for new row for large matrix
        row.append(h[0])                                        # Append mat1 id, as first item in each row
        for i in xrange(length-1):                              # For length of large matrix header (add 0 to each row) - header -1 for first '\t'
                row.extend('0')
        for n in xrange(1,49):                                  # Loop through each col id
                for k in mat2:                                  # For every row in mat2
                        if int(h[n]) == 1 and int(k[n]) == 1:   # If the pos (count for that particular col id) is 1 from mat1 and mat2 matrix;
                                pos = header.index(k[0])        # Get the position of the mat2 id
                                row[pos] = str(int(row[pos]) + 1)       # Add 1 to current position in row - [i][j] = [mat1_id][mat2_id]
        print >> out, '\t'.join(row)                            # When row is completed (All columns are compared from both mat1 and mat2 matrices; print final row to large matrix
        total += 1                                              # Update running total
        sys.stdout.write('\r\t' + str(total) + '/' + str(tvh))  # Print progress to screen
        sys.stdout.flush()

print '\n######################################################################################'
print 'Matrix complete.'
print '########################################################################################\n'

下面是 mat1 中 id 前 30 次迭代的概要分析：

######################################################################################
Generating matrix by counts from smaller matrices.
########################################################################################


Total mat1 rows: 108887
Total mat2 rows: 55482

Progress:
        30/108887^C         2140074 function calls in 101.234 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1   70.176   70.176  101.234  101.234 build_matrix.py:3(<module>)
        4    0.000    0.000    0.000    0.000 len
    55514    0.006    0.000    0.006    0.000 method 'append' of 'list' objects
        1    0.000    0.000    0.000    0.000 method 'disable' of '_lsprof.Profiler' objects
  1719942    1.056    0.000    1.056    0.000 method 'extend' of 'list' objects
       30    0.000    0.000    0.000    0.000 method 'flush' of 'file' objects
    35776   29.332    0.001   29.332    0.001 method 'index' of 'list' objects
       31    0.037    0.001    0.037    0.001 method 'join' of 'str' objects
   164370    0.589    0.000    0.589    0.000 method 'split' of 'str' objects
   164370    0.033    0.000    0.033    0.000 method 'strip' of 'str' objects
       30    0.000    0.000    0.000    0.000 method 'write' of 'file' objects
        2    0.000    0.000    0.000    0.000 next
        3    0.004    0.001    0.004    0.001 open

我还计算了每次迭代的时间，每个 mat1 id 大约需要 2.5-3 秒，如果我没记错的话，完成整个过程大约需要 90 小时。这是关于从头到尾运行整个脚本所需的内容。

我已经编辑了一些主要位，通过将“0”乘以标题的长度来删除通过附加和 xrange 生成行的步骤，从而一步生成列表。我还制作了一个 mat2 id 的字典，其中索引为值，认为 dict 查找会比索引更快..

headdict = 
for k,v in enumerate(header):
        headdict[v] = k

total = 0                                                       # Running total - for progress meter
for h in mat1:                                                  # Loop through all mat1 ids - each row in the HC matrix
        timestart = time.clock()
        row = [h[0]] + ['0']*(length-1)                 # Empty list for new row for large matrix
        #row.append(h[0])                                       # Append mat1 id, as first item in each row
        #for i in xrange(length-1):                             # For length of large matrix header (add 0 to each row) - header -1 for first '\t'
        #       row.append('0')
        for n in xrange(1,49):                                  # Loop through each col id
                for k in mat2:                                  # For every row in mat2
                        if int(h[n]) == 1 and int(k[n]) == 1:   # If the pos (count for that particular col id) is 1 from mat1 and mat2 matrix;
                                pos = headdict[k[0]] #header.index(k[0])        # Get the position of the mat2 id
                                row[pos] = str(int(row[pos]) + 1)       # Add 1 to current position in row - [i][j] = [mat1_id][mat2_id]
        print >> out, '\t'.join(row)                            # When row is completed (All columns are compared from both mat1 and mat2 matrices; print final row to large matrix
        total += 1                                              # Update running total
        sys.stdout.write('\r\t' + str(total) + '/' + str(totmat1))  # Print progress to screen
        #sys.stdout.flush()
        timeend = time.clock()
        print timestart - timeend

【问题讨论】：

虽然我看到您的代码效率低下，但没有什么特别突出的。有效优化它的唯一方法是首先对其进行概要分析，这是您无法做到的。幸运的是，这是一项相当简单的任务 - 请参阅 How can you profile a Python script? 结果会告诉您需要处理哪些部分才能获得最快的速度。

你的解释很不清楚。 “如果每个矩阵的 pos (col) 为 1”是什么意思？当您说“那么最终的矩阵计数将增加 1”时，您指的是哪个计数？该计数将存储在哪里？这听起来有点模糊，就像您只想将第一个矩阵乘以第二个矩阵的转置，但无法确定。

@user2357112- 我添加了一个我希望（快速）完成的小例子。每个矩阵的 col 的 pos 基于 id（xrange(1,49) 中的 n）。每个 mat1 id 的行最初由 0 填充。如果 mat2 中存在相同 col-id 的 1，则通过遍历所有 mat2 id（新矩阵中的标题），计数会增加 1。在给定 mat1 id 的情况下，这些值按行存储，并且在该 id 完成后，它被写入 outfile。对于每次迭代，我将所有 mat1 行按 id 与 col 标头的每个 mat2 id 进行比较。

@user3358205：最后一行应该是0 1 1 1而不是0 0 1 1吗？如果是这样，这看起来就像你只想要一个矩阵相乘。

你没有转置第二个矩阵。 （另外，获得 NumPy。那家伙的矩阵乘法可能比你的要好，但 NumPy 的胜过任何一个版本。）

【参考方案1】：

这只是一个矩阵乘法。您想将第一个矩阵乘以第二个矩阵的转置。如需高效的矩阵运算，请获取NumPy。

如果您将两个输入矩阵读入 dtype numpy.int8 的 NumPy 数组，那么计算很简单

m1.dot(m2.T)

这需要几分钟，顶。

【参考方案2】：

我不太明白这段代码的作用（单字母变量名没有帮助）。

我的建议：尽量减少在最内层循环中执行的操作次数。比如内层需要重新计算pos吗？

pos = header.index(k[0])

如果可以对嵌套循环 k、h 和 n 重新排序，则可以减少代价高昂的 list.index，这是一个 O(n) 操作。

【讨论】：

查看最新编辑以将列表索引更改为 dict 查找。我需要 k[0] 索引来获取标题中的位置，以了解在给定 mat1 id (h) 的情况下在大矩阵的行中更新哪个值。我需要为每个 mat1 id (h) 迭代所有 mat2 id (k)，在 mat1 和 mat2 中逐个 col 查找值相同的位置。如果它们是组合矩阵中的 [mat1][mat2] 位置，则值将增加 1。