numpy: как вычислить якобиану матрицу

Question

Я сделал пример точечного произведения двух матриц, и я хочу вычислить матрицу Якоби этого точечного произведения. Первая матрица имеет форму 4х3, а вторая матрица имеет форму 2х4. Получающаяся в результате матрица Якоби должна иметь форму (4x3x2x3), потому что я вычисляю ее по первой матрице.

Итак, допустим, что каждый элемент матрицы Якоби представляет собой DiDjAkAl, этот элемент будет представлять частную производную от i , вывод j относительно входа k, l. Вот пример в коде:

import numpy as np

A = np.array([[1,2,3,4],
              [5,6,7,8],
              [9,10,11,12]])
#shape = (4x3)

b = np.array([[1,2],
              [3,4],
              [5,6],
              [7,8]])
#shape = (2x4)

dotProduct = A.dot(b)
#shape = (2x3)

Так что в основном мой вопрос заключается в том, как найти якобиан этого точечного произведения.

Answer 1

Исходя из вашего решения, мы видим, что индексы i и k всегда являются одним и тем же значением, поэтому для l oop можно удалить. Тогда первое для l oop можно назначить всем в одном go, используя список в качестве слайса. Конкретно:

def vect_jacobian(A,b):

    dotProduct = A.dot(b)
    ans = np.zeros((A.shape[0], A.shape[1], dotProduct.shape[0], dotProduct.shape[1]))
    L = list(range(A.shape[0]))
    ans[L,:,L,:] = b
    return ans

r = vect_jacobian(A,b)

Используя %timeit на моем локальном компьютере, мы можем видеть, что даже с этими матрицами небольшого размера улучшение в 3 раза:

%timeit jacobian(A,b)
29.6 µs ± 3.74 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit vect_jacobian(A,b)
10.4 µs ± 589 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Answer 2

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

def jacobian(A, b):
    dotProduct = A.dot(b)
    ans = np.zeros((A.shape[0], A.shape[1], dotProduct.shape[0], dotProduct.shape[1]))
    for l in range(A.shape[0]):
        i = 0
        for k in range(A.shape[1]):
            ans[l][k][l][:] = b[i][:]
            i = i + 1
    return ans