计算 Python 中的雅可比矩阵
Compute the Jacobian matrix in 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 矩阵的内容，试图构建一个，根据我目前所读的内容，此 python 代码应被视为 jacobian。我理解的对吗？
雅可比行列式仅针对 vector-valued 函数 定义。您不能使用充满常量的数组来计算雅可比行列式；您必须知道基础函数及其偏导数，或它们的数值近似值。当您认为常数（相对于某物）的（偏）导数是 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
您可以使用哈佛 autograd 库 (link)，其中 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]))
否则，您可以使用可用于 sympy 中的矩阵的 jacobian 方法：
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)
此外，您可能也有兴趣查看此低级别变体 (link). MATLAB provides nice documentation on its jacobian function here。
更新：请注意，autograd 库已被纳入 jax，它提供了计算正向和逆雅可比矩阵的函数 (link)。
在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]])
这里是向量函数 f(x) 的数学雅可比行列式的 Python 实现，它被假定为 return 一维 numpy 数组。
import numpy as np
def J(f, x, dx=1e-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。
如果您想同时为多个点求出雅可比行列式的数值（例如，如果您的函数接受形状 (n, x) 并输出 (n, y)），这里有一个函数。这基本上是詹姆斯卡特的回答，但有很多要点。 dx 可能需要根据他的回答中的绝对值进行调整。
def numerical_jacobian(f, xs, dx=1e-6):
"""
f is a function that accepts input of shape (n_points, input_dim)
and outputs (n_points, output_dim)
return the jacobian as (n_points, output_dim, input_dim)
"""
if len(xs.shape) == 1:
xs = xs[np.newaxis, :]
assert len(xs.shape) == 2
ys = f(xs)
x_dim = xs.shape[1]
y_dim = ys.shape[1]
jac = np.empty((xs.shape[0], y_dim, x_dim))
for i in range(x_dim):
x_try = xs + dx * e(x_dim, i + 1)
jac[:, :, i] = (f(x_try) - ys) / dx
return jac
def e(n, i):
ret = np.zeros(n)
ret[i - 1] = 1.0
return ret
虽然 autograd 是一个很好的库，但请务必查看它的 upgraded version JAX，它有很好的文档记录（与 autograd 相比）。
一个简单的例子：
import jax.numpy as jnp
from jax import jacfwd
# Define some simple function.
def sigmoid(x):
return 0.5 * (jnp.tanh(x / 2) + 1)
# Note that here, I want a derivative of a "vector" output function (inputs*a + b is a vector) wrt a input
# "vector" a at a0: Derivative of vector wrt another vector is a matrix: The Jacobian
def simpleJ(a, b, inputs): #inputs is a matrix, a & b are vectors
return sigmoid(jnp.dot(inputs, a) + b)
inputs = jnp.array([[0.52, 1.12, 0.77],
[0.88, -1.08, 0.15],
[0.52, 0.06, -1.30],
[0.74, -2.49, 1.39]])
b = jnp.array([0.2, 0.1, 0.3, 0.2])
a0 = jnp.array([0.1,0.7,0.7])
# Isolate the function: variables to be differentiated from the constant parameters
f = lambda a: simpleJ(a, b, inputs) # Now f is just a function of variable to be differentiated
J = jacfwd(f)
# Till now I have only calculated the derivative, it still needs to be evaluated at a0.
J(a0)
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 矩阵的内容，试图构建一个，根据我目前所读的内容，此 python 代码应被视为 jacobian。我理解的对吗？
雅可比行列式仅针对 vector-valued 函数 定义。您不能使用充满常量的数组来计算雅可比行列式；您必须知道基础函数及其偏导数，或它们的数值近似值。当您认为常数（相对于某物）的（偏）导数是 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
您可以使用哈佛 autograd 库 (link)，其中 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]))
否则，您可以使用可用于 sympy 中的矩阵的 jacobian 方法：
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)
此外，您可能也有兴趣查看此低级别变体 (link). MATLAB provides nice documentation on its jacobian function here。
更新：请注意，autograd 库已被纳入 jax，它提供了计算正向和逆雅可比矩阵的函数 (link)。
在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]])
这里是向量函数 f(x) 的数学雅可比行列式的 Python 实现，它被假定为 return 一维 numpy 数组。
import numpy as np
def J(f, x, dx=1e-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。
如果您想同时为多个点求出雅可比行列式的数值（例如，如果您的函数接受形状 (n, x) 并输出 (n, y)），这里有一个函数。这基本上是詹姆斯卡特的回答，但有很多要点。 dx 可能需要根据他的回答中的绝对值进行调整。
def numerical_jacobian(f, xs, dx=1e-6):
"""
f is a function that accepts input of shape (n_points, input_dim)
and outputs (n_points, output_dim)
return the jacobian as (n_points, output_dim, input_dim)
"""
if len(xs.shape) == 1:
xs = xs[np.newaxis, :]
assert len(xs.shape) == 2
ys = f(xs)
x_dim = xs.shape[1]
y_dim = ys.shape[1]
jac = np.empty((xs.shape[0], y_dim, x_dim))
for i in range(x_dim):
x_try = xs + dx * e(x_dim, i + 1)
jac[:, :, i] = (f(x_try) - ys) / dx
return jac
def e(n, i):
ret = np.zeros(n)
ret[i - 1] = 1.0
return ret
虽然 autograd 是一个很好的库，但请务必查看它的 upgraded version JAX，它有很好的文档记录（与 autograd 相比）。
一个简单的例子：
import jax.numpy as jnp
from jax import jacfwd
# Define some simple function.
def sigmoid(x):
return 0.5 * (jnp.tanh(x / 2) + 1)
# Note that here, I want a derivative of a "vector" output function (inputs*a + b is a vector) wrt a input
# "vector" a at a0: Derivative of vector wrt another vector is a matrix: The Jacobian
def simpleJ(a, b, inputs): #inputs is a matrix, a & b are vectors
return sigmoid(jnp.dot(inputs, a) + b)
inputs = jnp.array([[0.52, 1.12, 0.77],
[0.88, -1.08, 0.15],
[0.52, 0.06, -1.30],
[0.74, -2.49, 1.39]])
b = jnp.array([0.2, 0.1, 0.3, 0.2])
a0 = jnp.array([0.1,0.7,0.7])
# Isolate the function: variables to be differentiated from the constant parameters
f = lambda a: simpleJ(a, b, inputs) # Now f is just a function of variable to be differentiated
J = jacfwd(f)
# Till now I have only calculated the derivative, it still needs to be evaluated at a0.
J(a0)