Эффективно вычислить коммутационную матрицу в numpy / scipy - PullRequest
Новая
       26

Эффективно вычислить коммутационную матрицу в numpy / scipy

1 голос
/ 14 марта 2020

Я пытаюсь вычислить матрицу коммутации в python для большого набора данных. Я написал следующий код, но обнаружил, что он работает ужасно (и встречается с ошибками памяти для примеров около 500 на 500). В моем коде a и b эквивалентны нотациям m и n на связанной странице википедии. Кто-нибудь может предложить более быструю и более эффективную память для моей текущей попытки? 

def vec(matrix):

    #Return vectorised matrix
    return(matrix.transpose().reshape(matrix.shape[0]*matrix.shape[1],1))


def commutation(a, b):

    # Example matrix with unique elements
    m = np.arange(a*b).reshape(a,b)

    # Vec(m) 
    vecm = vec(m)
    vecm = vecm.reshape(vecm.shape[0])

    # Get row inds
    rowInds = np.arange(a*b)

    # Get column inds
    colInds = np.argsort(vecm)
    colInds = colInds.reshape(colInds.shape[0])

    # Work out mapping between them.
    K = scipy.sparse.csr_matrix((np.ones(a*b),(rowInds,colInds)))

    return(K)

1 Ответ

2 голосов
/ 14 марта 2020

Ниже приведена улучшенная версия вашего кода: 

import numpy as np
from scipy.sparse import csr_matrix

def vec(A):
    m, n = A.shape[0], A.shape[1]
    return A.reshape(m*n, order='F')

def commutation_matrix_sp(A):
    m, n = A.shape[0], A.shape[1]
    row  = np.arange(m*n)
    col  = row.reshape((m, n), order='F').ravel()
    data = np.ones(m*n, dtype=np.int8)
    K = csr_matrix((data, (row, col)), shape=(m*n, m*n))
    return K

Тест: 

A = np.random.rand(500, 500)
K = commutation_matrix_sp(A)

print(f'{K.data.nbytes/2**20:.2f} MB')
# 0.24 MB

print(np.all(K @ vec(A) == vec(A.T)))
# True
...