Вычислить матрицу Якоби в Python

Вопрос:import numpy as np a = np.array([[1,2,3], [4,5,6], [7,8,9]]) b = np.array([[1,2,3]]).T c = a.dot(b) #function jacobian = a # as partial derivative of c w.r.t to b is a. Я читаю о jacobian Matrix, пытаясь построить один и из того, что я прочитал до сих пор, этот код python следует рассматривать как jacobian. Я

Вопрос:import numpy as np a = np.array([[1,2,3], [4,5,6], [7,8,9]]) b = np.array([[1,2,3]]).T c = a.dot(b) #function jacobian = a # as partial derivative of c w.r.t to b is a.

Я читаю о jacobian Matrix, пытаясь построить один и из того, что я прочитал до сих пор, этот код python следует рассматривать как jacobian. Я понимаю это право?

Ответ №1

Вы можете использовать autograd библиотеку Гарварда (ссылка), где grad и jacobian принимают функцию в качестве аргумента:

import autograd.numpy as np from autograd import grad, jacobian x = np.array([5,3], dtype=float) def cost(x): return x[0]**2 / x[1] — np.log(x[1]) gradient_cost = grad(cost) jacobian_cost = jacobian(cost) gradient_cost(x) jacobian_cost(np.array([x,x,x]))

В противном случае вы можете использовать метод jacobian доступный для матриц в sympy:

from sympy import sin, cos, Matrix from sympy.abc import rho, phi X = Matrix([rho*cos(phi), rho*sin(phi), rho**2]) Y = Matrix([rho, phi]) X.jacobian(Y)

Также вам может быть интересно посмотреть этот низкоуровневый вариант (ссылка).

Ответ №2

Якобиан определяется только для вектор- функций. Вы не можете работать с массивами, заполненными константами, для вычисления якобиана; вы должны знать основную функцию и ее частные производные или численное их приближение. Это очевидно, если учесть, что (частичная) производная от константы (по отношению к чему-то) равна 0.

В Python вы можете работать с символическими математическими модулями, такими как SymPy или SymEngine для вычисления якобиан функций. Здесь простая демонстрация примера из Википедии:

Использование модуля SymEngine:

Python 2.7.11 (v2.7.11:6d1b6a68f775, Dec 5 2015, 20:40:30) [MSC v.1500 64 bit (AMD64)] on win32 Type «help», «copyright», «credits» or «license» for more information. >>> >>> import symengine >>> >>> >>> vars = symengine.symbols(‘x y’) # Define x and y variables >>> f = symengine.sympify([‘y*x**2’, ‘5*x + sin(y)’]) # Define function >>> J = symengine.zeros(len(f),len(vars)) # Initialise Jacobian matrix >>> >>> # Fill Jacobian matrix with entries … for i, fi in enumerate(f): … for j, s in enumerate(vars): … J[i,j] = symengine.diff(fi, s) … >>> print J [2*x*y, x**2] [5, cos(y)] >>> >>> print symengine.Matrix.det(J) 2*x*y*cos(y) — 5*x**2 Ответ №3

В python 3 вы можете попробовать пакет sympy:

import sympy as sym def Jacobian(v_str, f_list): vars = sym.symbols(v_str) f = sym.sympify(f_list) J = sym.zeros(len(f),len(vars)) for i, fi in enumerate(f): for j, s in enumerate(vars): J[i,j] = sym.diff(fi, s) return J Jacobian(‘u1 u2’, [‘2*u1 + 3*u2′,’2*u1 — 3*u2’])

который выдает:

Matrix([[2, 3],[2, -3]]) Ответ №4

Вот Python-реализация математического якобиана вектор-функции f(x), которая, как предполагается, возвращает одномерный массив с нулевыми значениями.

import numpy as np def J(f, x, dx=10^-8): n = len(x) func = f(x) jac = np.zeros((n, n)) for j in range(n): #through columns to allow for vector addition Dxj = (abs(x[j])*dx if x[j] != 0 else dx) x_plus = [(xi if k != j else xi+Dxj) for k, xi in enumerate(x)] jac[:, j] = (f(x_plus)-func)/Dxj return jac

Рекомендуется сделать dx ~ 10 ^ -8.

Оцените статью
Добавить комментарий