Làm thế nào để Matlab transpose một ma trận thưa thớt?

Question

Tôi đã tự hỏi về câu hỏi này trong một thời gian nhưng không thể tìm thấy một tham chiếu: Làm thế nào để Matlab transpose một ma trận thưa thớt quá nhanh, cho rằng nó được lưu trữ trong CSC (nén thưa thớt cột) định dạng?

Cũng its documentation xác minh tính hiệu quả của chuyển vị ma trận thưa thớt:

Để làm điều này (truy cập vào từng hàng), bạn có thể transpose ma trận, thực hiện các hoạt động trên các cột, và sau đó retranspose kết quả ... Hiện cần thiết để chuyển đổi ma trận là không đáng kể.

Follow-up (sửa đổi theo đề nghị của @Mikhail):

Tôi đồng ý với @Roger và @Milhail rằng cài đặt một lá cờ là đủ cho nhiều hoạt động như các hoạt động BLAS hoặc BLAS thưa thớt về giao diện của họ. Nhưng nó xuất hiện với tôi rằng Matlab thực sự chuyển đổi "thực tế". Ví dụ, tôi có một ma trận X thưa thớt với kích thước m * n = 7984 * 12.411, và tôi muốn mở rộng mỗi cột và mỗi hàng:

% scaling each column 
t = 0; 
for i = 1 : 1000 
    A = X; t0 = tic; 
    A = bsxfun(@times, A, rand(1,n)); 
    t = t + toc(t0); 
end 

t = 0,023636 giây

% scaling each row 
t = 0; 
for i = 1 : 1000 
    A = X; t0 = tic; 
    A = bsxfun(@times, A, rand(m,1)); 
    t = t + toc(t0); 
end 

t = 138.3586 giây

% scaling each row by transposing X and transforming back 
t = 0; 
for i = 1 : 1000 
    A = X; t0 = tic; 
    A = A'; A = bsxfun(@times, A, rand(1,m)); A = A'; 
    t = t + toc(t0); 
end 

t = 19,5433 giây

kết quả này có nghĩa là truy cập vào cột theo cột nhanh hơn truy cập hàng theo hàng. Nó có ý nghĩa bởi vì ma trận thưa thớt được lưu trữ theo từng cột. Vì vậy, lý do duy nhất cho tốc độ nhanh của tỷ lệ cột của X 'nên X thực sự được chuyển sang X' thay vì đặt cờ.

Ngoài ra, nếu mọi ma trận thưa thớt được lưu trữ ở định dạng CSC, chỉ cần đặt cờ không thể tạo X 'ở định dạng CSC.

Mọi nhận xét? Cảm ơn trước.

Answer 1

Tôi đồng ý với những gì Roger Rowland đã đề cập trong nhận xét. Để đưa ra đề xuất này, bạn có thể kiểm tra một số chức năng từ giao diện BLAS, mà MATLAB sử dụng cho các phép toán ma trận. Tôi không chắc nó sử dụng công cụ gì, nhưng vì họ sử dụng Intel IPP để xử lý hình ảnh, tôi cho rằng họ cũng có thể sử dụng Intel MKL để làm cho các hoạt động ma trận nhanh chóng.

Và đây là tài liệu cho hàm mkl_?cscsv, giải quyết một hệ phương trình tuyến tính cho ma trận thưa thớt theo định dạng CSC. Lưu ý cờ đầu vào transa, trong đó xác định rõ liệu ma trận được cung cấp có được coi là được hoán vị hay không.

Answer 2

Sau khi khám phá một tuần, tôi đoán về cơ chế nội bộ của việc chuyển đổi ma trận đang sắp xếp.

Giả sử A là một ma trận thưa thớt,

[I, J, S] = find(A); 
[sorted_I, idx] = sort(I); 
J = J(idx); 
S = S(idx); 
B = sparse(J, sorted_I, S); 

Sau đó B là transpose của A.

Việc triển khai ở trên có khoảng một nửa hiệu quả của tích hợp sẵn của máy Matlab transpose trên máy của tôi.Xem xét các chức năng tích hợp của Matlab là đa luồng, phỏng đoán của tôi có thể là hợp lý.

Answer 3

Tôi nhận ra mình hơi muộn trong trò chơi, nhưng tôi nghĩ tôi có thể giúp giải thích một chút về câu hỏi này. Transposing một ma trận thưa thớt thực sự là một nhiệm vụ đơn giản có thể được thực hiện trong thời gian tỷ lệ thuận với số lượng các phần tử nonzero trong ma trận đầu vào. Giả sử rằng A là một ma trận mxn được lưu trữ ở định dạng CSC, tức là, A được xác định bởi ba mảng:

1. elemsA, chiều dài NNZ (A), mà các cửa hàng các yếu tố khác không A
2. prowA, chiều dài nnz (A), lưu trữ các chỉ số hàng của các phần tử không đồng nhất trong A
3. pcolA, có độ dài n + 1, sao cho tất cả các phần tử không phải trong cột j của A được lập chỉ mục theo phạm vi [pcolA (j), pcolA (j + 1))

Nếu B biểu thị chuyển vị A, thì mục tiêu của chúng tôi là xác định các mảng tương tự elemsB, prowB, pcolB. Để làm như vậy, chúng tôi sử dụng thực tế là các hàng của A tạo thành các cột của B. Let tmp là một mảng sao cho tmp (1) = 0 và tmp (i + 1) là số phần tử trong hàng i của A cho i = 1, ..., m. Sau đó, nó theo sau đó tmp (i + 1) là số lượng các phần tử trong cột i của B. Do đó, tổng tích lũy của tmp là giống như pcolB. Bây giờ giả sử rằng tmp đã được ghi đè bằng tổng tích lũy của nó. Sau đó elemsB và prowB có thể được dân cư như sau

for j = 1,...,n 
     for k = pcolA(j),...,pcolA(j + 1) - 1 
      prowB(tmp(prowA(k))) = j 
      elemsB(tmp(prowA(k))) = elemsA(k) 
      tmp(prowA(k)) = tmp(prowA(k)) + 1 
     end 
    end 

Các tmp mảng được sử dụng để chỉ mục vào prowB và elemsB khi thêm một yếu tố mới và sau đó được cập nhật cho phù hợp. Đưa này hoàn toàn, chúng ta có thể viết một file mex trong C++ mà thực hiện các thuật toán chuyển vị:

#include "mex.h" 
#include <vector> 
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {  
    // check input output 
    if (nrhs != 1) 
     mexErrMsgTxt("One input argument required"); 
    if (nlhs > 1) 
     mexErrMsgTxt("Too many output arguments"); 

    // get input sparse matrix A 
    if (mxIsEmpty(prhs[0])) { // is A empty? 
     plhs[0] = mxCreateSparse(0, 0, 0, mxREAL); 
     return; 
    } 
    if (!mxIsSparse(prhs[0]) || mxIsComplex(prhs[0])) // is A real and sparse? 
     mexErrMsgTxt("Input matrix must be real and sparse"); 
    double* A = mxGetPr(prhs[0]);   // real vector for A 
    mwIndex* prowA = mxGetIr(prhs[0]);  // row indices for elements of A 
    mwIndex* pcolindexA = mxGetJc(prhs[0]); // index into the columns 
    mwSize M = mxGetM(prhs[0]);    // number of rows in A 
    mwSize N = mxGetN(prhs[0]);    // number of columns in A 

    // allocate memory for A^T 
    plhs[0] = mxCreateSparse(N, M, pcolindexA[N], mxREAL); 
    double* outAt = mxGetPr(plhs[0]); 
    mwIndex* outprowAt = mxGetIr(plhs[0]); 
    mwIndex* outpcolindexAt = mxGetJc(plhs[0]); 

    // temp[j + 1] stores the number of nonzero elements in row j of A 
    std::vector<mwSize> temp(M + 1, 0); 
    for(mwIndex i = 0; i != N; ++i) { 
     for(mwIndex j = pcolindexA[i]; j < pcolindexA[i + 1]; ++j) 
      ++temp[prowA[j] + 1]; 
    } 
    outpcolindexAt[0] = 0; 
    for(mwIndex i = 1; i <= M; ++i) { 
     outpcolindexAt[i] = outpcolindexAt[i - 1] + temp[i]; 
     temp[i] = outpcolindexAt[i]; 
    } 
    for(mwIndex i = 0; i != N; ++i) { 
     for(mwIndex j = pcolindexA[i]; j < pcolindexA[i + 1]; ++j) { 
      outprowAt[temp[prowA[j]]] = i; 
      outAt[temp[prowA[j]]++] = A[j]; 
     } 
    } 
} 

So sánh thuật toán này với việc thực hiện các chuyển vị Matlab, chúng tôi quan sát thời gian thực hiện tương tự. Lưu ý rằng thuật toán này có thể được sửa đổi một cách đơn giản để loại bỏ mảng temp.