Tìm ma trận tương quan

Question

Tôi có ma trận khá lớn (khoảng 50K hàng) và tôi muốn in hệ số tương quan giữa mỗi hàng trong ma trận. Tôi đã viết mã Python như thế này:

for i in xrange(rows): # rows are the number of rows in the matrix. 
    for j in xrange(i, rows): 
     r = scipy.stats.pearsonr(data[i,:], data[j,:]) 
     print r 

Xin lưu ý rằng tôi đang làm cho việc sử dụng pearsonr chức năng có sẵn từ các mô-đun scipy (http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html).

Câu hỏi của tôi là: Có cách nào nhanh hơn để thực hiện việc này không? Có một số kỹ thuật phân vùng ma trận mà tôi có thể sử dụng không?

Cảm ơn!

Answer 1

bạn có thể sử dụng các mô-đun python đa tiến, đoạn lên hàng của bạn thành 10 bộ, đệm kết quả của bạn và sau đó in những thứ ra (điều này sẽ chỉ đẩy nó lên trên một máy đa lõi dù)

http://docs.python.org/library/multiprocessing.html

btw: bạn cũng phải chuyển đoạn mã của mình thành một hàm và cũng xem xét cách thực hiện khôi phục dữ liệu. có mỗi tiến trình con có một danh sách như thế này ... [startcord, stopcord, buff] .. có thể làm việc độc đáo

def myfunc(thelist): 
    for i in xrange(thelist[0]:thelist[1]): 
    .... 
    thelist[2] = result 

Answer 2

New Solution

Sau khi xem xét câu trả lời của Joe Kington, tôi quyết định để xem xét mã corrcoef() và được lấy cảm hứng từ mã đó để thực hiện việc triển khai sau.

ms = data.mean(axis=1)[(slice(None,None,None),None)] 
datam = data - ms 
datass = np.sqrt(scipy.stats.ss(datam,axis=1)) 
for i in xrange(rows): 
    temp = np.dot(datam[i:],datam[i].T) 
    rs = temp/(datass[i:]*datass[i]) 

Mỗi vòng lặp tạo ra hệ số Pearson giữa hàng i và hàng i đến hàng cuối cùng. Nó rất nhanh. Đó là ít nhất 1.5x nhanh như việc sử dụng corrcoef() một mình vì nó không tính toán thừa các hệ số và một vài thứ khác. Nó cũng sẽ nhanh hơn và sẽ không cung cấp cho bạn các vấn đề về bộ nhớ với ma trận hàng 50.000 vì sau đó bạn có thể chọn lưu trữ từng bộ r hoặc xử lý chúng trước khi tạo một bộ khác. Nếu không lưu trữ bất kỳ dài hạn của r, tôi đã có thể nhận được mã trên để chạy trên 50.000 x 10 bộ dữ liệu được tạo ngẫu nhiên trong một phút trên máy tính xách tay khá mới của tôi.

Cũ Solution

Trước tiên, tôi sẽ không khuyên bạn nên in ra r của màn hình. Đối với 100 hàng (10 cột), đây là sự khác biệt 19,79 giây với in so với 0,01 giây mà không sử dụng mã của bạn. Chỉ cần lưu trữ của r và sử dụng chúng sau này nếu bạn muốn, hoặc làm một số chế biến trên chúng khi bạn đi cùng như tìm kiếm một số r lớn nhất.

Thứ hai, bạn có thể nhận được một số khoản tiết kiệm bằng cách không tính thừa một số lượng. Hệ số Pearson được tính toán trong scipy bằng cách sử dụng một số đại lượng mà bạn có thể tính toán trước thay vì tính toán mỗi khi một hàng được sử dụng. Ngoài ra, bạn không sử dụng các giá trị p (mà cũng được trả về bởi pearsonr() vì vậy hãy làm xước rằng quá Sử dụng mã bên dưới:.

r = np.zeros((rows,rows)) 
ms = data.mean(axis=1) 

datam = np.zeros_like(data) 
for i in xrange(rows): 
    datam[i] = data[i] - ms[i] 
datass = scipy.stats.ss(datam,axis=1) 
for i in xrange(rows): 
    for j in xrange(i,rows): 
     r_num = np.add.reduce(datam[i]*datam[j]) 
     r_den = np.sqrt(datass[i]*datass[j]) 
     r[i,j] = min((r_num/r_den), 1.0) 

tôi nhận được một tốc độ lên khoảng 4.8x so với scipy thẳng mã khi tôi đã loại bỏ các công cụ giá trị p - 8.8x nếu tôi để lại các công cụ giá trị p trong đó (tôi đã sử dụng 10 cột với hàng trăm hàng). Tôi cũng kiểm tra xem nó có đưa ra kết quả tương tự không. một cải tiến thực sự rất lớn, nhưng nó có thể hữu ích.

Cuối cùng, bạn đang mắc kẹt với vấn đề mà bạn đang tính toán (50000) * (50001)/2 = 1,250,025.000 hệ số Pearson (nếu tôi đếm chính xác). Đó là rất nhiều. Bằng cách này, thực sự không cần tính toán hệ số Pearson của mỗi hàng với chính nó (nó sẽ bằng 1), nhưng điều đó chỉ giúp bạn tiết kiệm từ việc tính toán 50.000 hệ số Pearson. Với đoạn mã trên, tôi hy vọng rằng sẽ mất khoảng 4 1/4 giờ để tính toán nếu bạn có 10 cột dữ liệu dựa trên kết quả của tôi trên các tập dữ liệu nhỏ hơn.

Bạn có thể nhận được một số cải tiến bằng cách lấy mã trên vào Cython hoặc một cái gì đó tương tự. Tôi hy vọng rằng bạn sẽ có thể nhận được lên đến một cải tiến 10x trên thẳng Scipy nếu bạn may mắn. Ngoài ra, theo đề xuất của pyInTheSky, bạn có thể thực hiện một số xử lý đa.

Answer 3

Bạn đã thử chỉ sử dụng numpy.corrcoef? Thấy như bạn không sử dụng giá trị p, nó nên làm chính xác những gì bạn muốn, với càng ít phiền toái càng tốt. (Trừ khi tôi nhớ nhầm chính xác R của pearson là gì, điều này hoàn toàn có thể.)

Chỉ cần nhanh chóng kiểm tra kết quả trên dữ liệu ngẫu nhiên, nó trả về chính xác mã giống như @Justin Peel ở trên và chạy nhanh hơn 100 lần .

Ví dụ, thử nghiệm mọi thứ với 1000 hàng và 10 cột dữ liệu ngẫu nhiên ...:

import numpy as np 
import scipy as sp 
import scipy.stats 

def main(): 
    data = np.random.random((1000, 10)) 
    x = corrcoef_test(data) 
    y = justin_peel_test(data) 
    print 'Maximum difference between the two results:', np.abs((x-y)).max() 
    return data 

def corrcoef_test(data): 
    """Just using numpy's built-in function""" 
    return np.corrcoef(data) 

def justin_peel_test(data): 
    """Justin Peel's suggestion above""" 
    rows = data.shape[0] 

    r = np.zeros((rows,rows)) 
    ms = data.mean(axis=1) 

    datam = np.zeros_like(data) 
    for i in xrange(rows): 
     datam[i] = data[i] - ms[i] 
    datass = sp.stats.ss(datam,axis=1) 
    for i in xrange(rows): 
     for j in xrange(i,rows): 
      r_num = np.add.reduce(datam[i]*datam[j]) 
      r_den = np.sqrt(datass[i]*datass[j]) 
      r[i,j] = min((r_num/r_den), 1.0) 
      r[j,i] = r[i,j] 
    return r 

data = main() 

sản lượng một sự khác biệt tuyệt đối tối đa ~ 3.3e-16 giữa hai kết quả

Và timings :

In [44]: %timeit corrcoef_test(data) 
10 loops, best of 3: 71.7 ms per loop 

In [45]: %timeit justin_peel_test(data) 
1 loops, best of 3: 6.5 s per loop 

numpy.corrcoef chỉ nên làm những gì bạn muốn và nhanh hơn rất nhiều.