Матричные операции

Question

Я хочу вычислить следующие значения для всех i и j:

M_ki = Sum[A_ij - A_ik - A_kj + A_kk, 1 <= j <= n]

Как я могу сделать это, используя Numpy (Python) без явного цикла?

Спасибо!

Answer 1

Вот общая стратегия решения такого рода проблем.

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

import numpy as np
from numpy import newaxis

def explicit(a):
    n = a.shape[0]
    m = np.zeros_like(a)
    for k in range(n):
        for i in range(n):
            for j in range(n):
                m[k,i] += a[i,j] - a[i,k] - a[k,j] + a[k,k]
    return m

def implicit(a):
    n = a.shape[0]
    m = np.zeros_like(a)
    for k in range(n):
        for i in range(n):
            for j in range(n):
                m[k,i] += a[i,j] - a[i,k] - a[k,j] + a[k,k]
    return m

a = np.random.randn(10,10)
assert np.allclose(explicit(a), implicit(a), atol=1e-10, rtol=0.)

Затем шаг за шагом векторизируйте функцию, отредактировав implicit, запуская скрипт на каждом шаге, чтобы убедиться, что они продолжают оставаться прежними:

Шаг 1

def implicit(a):
    n = a.shape[0]
    m = np.zeros_like(a)
    for k in range(n):
        for i in range(n):
            m[k,i] = (a[i,:] - a[k,:]).sum() - n*a[i,k] + n*a[k,k]
    return m

Шаг 2

def implicit(a):
    n = a.shape[0]
    m = np.zeros_like(a)
    m = - n*a.T + n*np.diag(a)[:,newaxis]
    for k in range(n):
        for i in range(n):
            m[k,i] += (a[i,:] - a[k,:]).sum()
    return m

Шаг 3

def implicit(a):
    n = a.shape[0]
    m = np.zeros_like(a)
    m = - n*a.T + n*np.diag(a)[:,newaxis]
    m += (a.T[newaxis,...] - a[...,newaxis]).sum(1)
    return m

Et voila '! Никаких петель в последнем. Чтобы векторизовать такие уравнения, вещание - это путь!

Предупреждение: убедитесь, что explicit - это уравнение, которое вы хотели векторизовать. Я не был уверен, должны ли суммы, которые не зависят от j, также суммироваться.