numpy более эффективная подматрица

Question

В настоящее время я использую приведенный ниже код, чтобы удалить подматрицу с i-й строкой и j-м столбцом, но после профилирования моего кода это, по-видимому, является одним из основных узких мест в моем коде. Есть ли более эффективный способ?

def submatrix(A, i, j): 
    logger.debug('submatrix(%r, %r, %r)', A, i, j) 
    B = empty(shape=tuple(x - 1 for x in A.shape), dtype=int) 
    B[:i, :j] = A[:i, :j] 
    B[i:, :j] = A[i+1:, :j] 
    B[:i, j:] = A[:i, j+1:] 
    B[i:, j:] = A[i+1:, j+1:] 
    return B 

     25015049 function calls (24599369 primitive calls) in 44.587 seconds 

    Ordered by: internal time 

    ncalls tottime percall cumtime percall filename:lineno(function) 
    3983040 15.541 0.000 20.719 0.000 defmatrix.py:301(__getitem__) 
    415680 10.216 0.000 33.069 0.000 hill.py:127(submatrix) 
415686/6 3.232 0.000 44.578 7.430 hill.py:112(det) 

Edit: Jaime предоставил хороший способ аппроксимировать модульные инверсий с использованием регулярных обратной и детерминанта, однако с большими базами (по модулю 256 в моем случае), неточность достаточно, чтобы оказать всю вещь спорным. Основное время раковина оказывается на самом деле быть GetItem в NumPy, но я считаю, что это вызван этими линиями:

B[:i, :j] = A[:i, :j] 
    B[i:, :j] = A[i+1:, :j] 
    B[:i, j:] = A[:i, j+1:] 
    B[i:, j:] = A[i+1:, j+1:] 

Это возможно бутылочное горлышко не копируя матрицы вокруг в памяти, но доступ матричного элемента ,

Answer 1

Хм ... вы только копируете матрицы, поэтому, вероятно, будет сложно ускорить работу, но одна вещь, которую вы можете попробовать, - проверить, что A находится в непрерывном блоке памяти, что может ускорить доступ к C. Посмотрите на numpy.ascontiguousarray().

Answer 2

Насколько я могу судить, submatrix просто удаляет i-й ряд и j-й столбец. Вы можете сделать это с np.delete

i = 3 
j = 4 
a = np.arange(100).reshape(10,10) 
b = np.delete(np.delete(a, i, 0), j, 1) 

Но для Расон @Jaime ссылается, это на самом деле медленнее :-/

timeit submatrix(a, i, j) 
#10000 loops, best of 3: 23.2 us per loop 

timeit subdel(a, i, j) 
#10000 loops, best of 3: 42.6 us per loop 

Но я оставлю это здесь сейчас.

Answer 3

Вычисление инверсии матрицы с использованием детерминант - очень медленный подход, независимо от того, как вы выполняете подматрицы. Давайте принимать глупый пример:

a = np.array([[3, 0, 2], 
       [2, 0, -2], 
       [0, 1, 1]]) 

Вы можете вычислить обратную быстро, как:

>>> np.linalg.inv(a) 
array([[ 0.2, 0.2, 0. ], 
     [-0.2, 0.3, 1. ], 
     [ 0.2, -0.3, -0. ]]) 

Но вычислить модульную обратный, вам необходимо иметь как целочисленную матрицу, разделенную на целое число раза. Это целое число фактор, конечно, будет определяющим, так что вы можете сделать следующее:

>>> np.linalg.inv(a) * np.linalg.det(a) 
array([[ 2., 2., 0.], 
     [ -2., 3., 10.], 
     [ 2., -3., -0.]]) 

И обратное a это целая матрица, делится на определитель a. В качестве функции, вы можете сделать:

def extended_euclidean_algorithm(a, b) : 
    """ 
    Computes a solution to a x + b y = gcd(a,b), as well as gcd(a,b), 
    using the extended Euclidean algorithm. 
    """ 
    if b == 0 : 
     return 1, 0, a 
    else : 
     x, y, gcd = extended_euclidean_algorithm(b, a % b) 
     return y, x - y * (a // b), gcd 

def mmi(a, m) : 
    """ 
    Computes the modular multiplicative inverse of a modulo m, using the 
    extended Euclidean algorithm. 
    """ 
    x, y, gcd = extended_euclidean_algorithm(a, m) 
    if gcd == 1 : 
     return x % m 
    else : 
     return None 

def modular_inv(a, m): 
    det_a = np.linalg.det(a) 
    inv_a = np.linalg.inv(a) * det_a 
    det_a = np.rint(det_a).astype(int) 
    inv_a = np.rint(inv_a).astype(int) 
    return ((inv_a % m) * mmi(det_a, m)) % m 

А теперь:

>>> a = np.random.randint(10, size=(10, 10)) 
>>> b = modular_inv(a, 7) 
>>> a.dot(b) % 7 
array([[1, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 1, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 1, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 1, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 1, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 0, 1, 0, 0, 0, 0], 
     [0, 0, 0, 0, 0, 0, 1, 0, 0, 0], 
     [0, 0, 0, 0, 0, 0, 0, 1, 0, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 1, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 1]])