nhân ma trận numpy để lưu trữ hình tam giác/thưa thớt?

Question

Tôi đang làm việc với vấn đề nhân ma trận thưa thớt rất lớn (matmul). Ví dụ: giả sử:

• A là ma trận nhị phân (75 x 200.000). Đó là thưa thớt, vì vậy tôi đang sử dụng csc để lưu trữ. Tôi cần phải làm các hoạt động matmul sau:

• B = A.transpose() * Một

• Kết quả sẽ là một ma trận thưa thớt và đối xứng của kích thước 200Kx200K.

Thật không may, B sẽ là cách đến lớn để lưu trữ trong RAM (hoặc "trong lõi") trên máy tính xách tay của tôi. Mặt khác, tôi may mắn vì có một số thuộc tính cho B nên giải quyết vấn đề này.

Vì B sẽ đối xứng dọc theo đường chéo và thưa thớt, tôi có thể sử dụng ma trận tam giác (trên/dưới) để lưu trữ kết quả hoạt động matmul và định dạng lưu trữ ma trận thưa thớt có thể giảm kích thước.

Câu hỏi của tôi là ... có thể bị sần sùi hoặc scipy, trước tiên, yêu cầu lưu trữ đầu ra sẽ trông như thế nào để tôi có thể chọn giải pháp lưu trữ bằng cách sử dụng gọn gàng và tránh "ma trận quá lớn" lỗi thời gian chạy sau vài phút (giờ) tính toán?

Nói cách khác, các yêu cầu lưu trữ cho ma trận nhân được xấp xỉ bằng cách phân tích nội dung của hai ma trận đầu vào bằng cách sử dụng thuật toán đếm gần đúng?

• https://en.wikipedia.org/wiki/Approximate_counting_algorithm

Nếu không, tôi đang nhìn vào một giải pháp brute force. Một cái gì đó liên quan đến bản đồ/giảm, out-of-core lưu trữ, hoặc một giải pháp phân matmul (strassens thuật toán) từ các liên kết web sau đây:

Một vài Map/Giảm các giải pháp chia nhỏ vấn đề

• http://www.norstad.org/matrix-multiply/index.html
• http://bpgergo.blogspot.com/2011/08/matrix-multiplication-in-python.html

Một giải pháp (PyTables) lưu trữ out-of-core

• Very large matrices using Python and NumPy

Một giải pháp phân matmul:

• https://en.wikipedia.org/wiki/Strassen_algorithm
• http://facultyfp.salisbury.edu/taanastasio/COSC490/Fall03/Lectures/FoxMM/example.pdf
• http://eli.thegreenplace.net/2012/01/16/python-parallelizing-cpu-bound-tasks-with-multiprocessing/

Cảm ơn trước cho bất cứ đề nghị, ý kiến, hoặc hướng dẫn!

Answer 1

Vì bạn sau khi sản phẩm của ma trận với vị trí của nó, giá trị tại [m, n] về cơ bản sẽ là sản phẩm chấm của các cột m và n trong ma trận ban đầu của bạn.

tôi sẽ sử dụng các ma trận sau đây là một ví dụ đồ chơi

a = np.array([[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1], 
       [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0], 
       [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1]]) 
>>> np.dot(a.T, a) 
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0], 
     [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1], 
     [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, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 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], 
     [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0], 
     [0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 2]]) 

Nó có hình dạng (3, 12) và có 7 mục khác không. Sản phẩm của nó transpose với nó là tất nhiên của hình dạng (12, 12) và có 16 số không mục, 6 của nó trong đường chéo, vì vậy nó chỉ đòi hỏi lưu trữ 11 yếu tố.

Bạn có thể có được một ý tưởng tốt về những gì kích thước của ma trận đầu ra của bạn sẽ là một trong hai cách sau:

CSR FORMAT

Nếu ma trận ban đầu của bạn có C cột khác không , ma trận mới của bạn sẽ có tối đa C**2 mục khác 0, trong đó C nằm trong đường chéo và được đảm bảo không bằng 0 và các mục còn lại bạn chỉ cần giữ một nửa, sao cho tối đa là (C**2 + C)/2 không yếu tố không. Tất nhiên, nhiều người trong số này cũng sẽ bằng không, vì vậy đây có lẽ là một đánh giá quá tổng.

Nếu ma trận của bạn được lưu trữ ở định dạng csr thì indices thuộc tính của scipy đối tượng tương ứng có một mảng với chỉ số cột của tất cả các zero yếu tố phi, vì vậy bạn có thể dễ dàng tính toán các ước tính trên như:

>>> a_csr = scipy.sparse.csr_matrix(a) 
>>> a_csr.indices 
array([ 2, 11, 1, 7, 10, 4, 11]) 
>>> np.unique(a_csr.indices).shape[0] 
6 

vì vậy, có 6 cột với mục khác không, và do đó ước tính sẽ là ít nhất 36 mục khác không, cách hơn so với thực 16.

CSC FORMAT

Nếu thay vì chỉ mục cột của các phần tử khác 0, chúng tôi có chỉ số hàng, chúng tôi thực sự có thể ước tính tốt hơn. Đối với sản phẩm chấm của hai cột không khác, chúng phải có phần tử khác 0 trong cùng một hàng. Nếu có R các phần tử khác 0 trong một hàng nhất định, chúng sẽ đóng góp R**2 các phần tử khác 0 vào sản phẩm. Khi bạn tính tổng số này cho tất cả các hàng, bạn buộc phải đếm một số phần tử nhiều lần, vì vậy đây cũng là một giới hạn trên.

Các chỉ số hàng của các yếu tố khác không của ma trận của bạn đang trong indices thuộc tính của một ma trận csc thưa thớt, vì vậy ước tính này có thể được tính như sau:

>>> a_csc = scipy.sparse.csc_matrix(a) 
>>> a_csc.indices 
array([1, 0, 2, 1, 1, 0, 2]) 
>>> rows, where = np.unique(a_csc.indices, return_inverse=True) 
>>> where = np.bincount(where) 
>>> rows 
array([0, 1, 2]) 
>>> where 
array([2, 3, 2]) 
>>> np.sum(where**2) 
17 

Đây là darn gần với thực 16! Và nó thực sự không phải là một sự trùng hợp rằng ước tính này thực sự là giống như:

>>> np.sum(np.dot(a.T,a),axis=None) 
17 

Trong mọi trường hợp, các mã sau đây sẽ cho phép bạn để thấy rằng việc ước tính là khá tốt:

def estimate(a) : 
    a_csc = scipy.sparse.csc_matrix(a) 
    _, where = np.unique(a_csc.indices, return_inverse=True) 
    where = np.bincount(where) 
    return np.sum(where**2) 

def test(shape=(10,1000), count=100) : 
    a = np.zeros(np.prod(shape), dtype=int) 
    a[np.random.randint(np.prod(shape), size=count)] = 1 
    print 'a non-zero = {0}'.format(np.sum(a)) 
    a = a.reshape(shape) 
    print 'a.T * a non-zero = {0}'.format(np.flatnonzero(np.dot(a.T, 
                   a)).shape[0]) 
    print 'csc estimate = {0}'.format(estimate(a)) 

>>> test(count=100) 
a non-zero = 100 
a.T * a non-zero = 1065 
csc estimate = 1072 
>>> test(count=200) 
a non-zero = 199 
a.T * a non-zero = 4056 
csc estimate = 4079 
>>> test(count=50) 
a non-zero = 50 
a.T * a non-zero = 293 
csc estimate = 294