{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 使用函数变换计算雅可比矩阵和黑塞矩阵\n", "\n", "[![下载Notebook](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/master/resource/_static/logo_notebook.svg)](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/notebook/master/tutorials/experts/zh_cn/func_programming/mindspore_Jacobians_Hessians.ipynb) [![下载样例代码](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/master/resource/_static/logo_download_code.svg)](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/notebook/master/tutorials/experts/zh_cn/func_programming/mindspore_Jacobians_Hessians.py) [![查看源文件](https://mindspore-website.obs.cn-north-4.myhuaweicloud.com/website-images/master/resource/_static/logo_source.svg)](https://gitee.com/mindspore/docs/blob/master/tutorials/experts/source_zh_cn/func_programming/Jacobians_Hessians.ipynb)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 雅可比矩阵\n", "\n", "在介绍MindSpore提供的计算雅可比矩阵的方法之前，首先对雅可比矩阵进行介绍。\n", "\n", "我们首先定义一个映射$\\textbf{f}$:\n", "\n", "$$R^{n} \\longrightarrow R^{m}$$\n", "\n", "我们在这里使用符号$\\longmapsto$表示集合元素之间的映射，并加粗所有代表向量的符号，因此有：\n", "\n", "$$\\textbf{x} \\longmapsto \\textbf{f}(\\textbf{x})$$\n", "\n", "其中$\\textbf{x} = (x_{1}, x_{2},\\dots, x_{n})$，$\\textbf{f(x)} = (f_{1}(\\textbf{x}), f_{2}(\\textbf{x}),\\dots,f_{m}(\\textbf{x}))$。\n", "\n", "我们将所有从$R^{n}$到$R^{m}$的映射构成的集合记为$F_{n}^{m}$。\n", "在这里，我们将从一个函数（映射）集到另一个函数（映射）集的映射称为一个操作。易得，梯度操作$\\nabla$是在函数集$F_{n}^{1}$上同时进行n次偏导数的操作，将其定义为一个从函数集$F_{n}^{1}$到$F_{n}^{n}$的映射：\n", "\n", "$$\\nabla：F_{n}^{1} \\longrightarrow F_{n}^{n}$$\n", "\n", "广义的梯度操作$\\partial$，被定义为在函数集$F_{n}^{m}$上同时进行n次偏导数的操作，\n", "\n", "$$\\partial: F_{n}^{m} \\longrightarrow F_{n}^{m \\times n}$$\n", "\n", "雅可比矩阵就是将操作$\\partial$作用于$\\textbf{f}$后得到的结果，即\n", "\n", "$$\\textbf{f} \\longmapsto \\partial \\textbf{f} = (\\frac{\\partial \\textbf{f}}{\\partial x_{1}}, \\frac{\\partial \\textbf{f}}{\\partial x_{2}}, \\dots, \\frac{\\partial \\textbf{f}}{\\partial x_{n}})$$\n", "\n", "得到雅可比矩阵，\n", "\n", "$$J_{f} = \\begin{bmatrix}\n", " \\frac{\\partial f_{1}}{\\partial x_{1}} &\\frac{\\partial f_{1}}{\\partial x_{2}} &\\dots &\\frac{\\partial f_{1}}{\\partial x_{n}} \\\\\n", " \\frac{\\partial f_{2}}{\\partial x_{1}} &\\frac{\\partial f_{2}}{\\partial x_{2}} &\\dots &\\frac{\\partial f_{2}}{\\partial x_{n}} \\\\\n", " \\vdots &\\vdots &\\ddots &\\vdots \\\\\n", " \\frac{\\partial f_{m}}{\\partial x_{1}} &\\frac{\\partial f_{m}}{\\partial x_{2}} &\\dots &\\frac{\\partial f_{m}}{\\partial x_{n}}\n", "\\end{bmatrix}$$\n", "\n", "雅可比矩阵的应用：在自动微分的前向模式中，每一次前向传播，可以求出雅可比矩阵的一列。在自动微分的反向模式中，每一次反向传播，我们可以计算出雅可比矩阵的一行。\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 计算雅可比矩阵\n", "\n", "使用标准的自动微分系统很难高效地计算雅可比矩阵，但MindSpore却提供了能够高效计算雅可比的方法，下面对这些方法进行介绍。\n", "\n", "首先，我们定义一个函数`forecast`，该函数是一个简单的线性函数，并带有一个非线性激活函数。" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import time\n", "import mindspore\n", "from mindspore import ops\n", "from mindspore import jacrev, jacfwd, vmap, vjp, jvp, grad\n", "import numpy as np\n", "\n", "mindspore.set_seed(1)\n", "\n", "def forecast(weight, bias, x):\n", " return ops.dense(x, weight, bias).tanh()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "接下来，我们构造一些数据：一个权重张量`weight`，一个偏差张量`bias`，还有一个输入向量`x`。" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "D = 16\n", "weight = ops.randn(D, D)\n", "bias = ops.randn(D)\n", "x = ops.randn(D)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "函数`forecast`对输入向量`x`做如下映射变换，$R^{D}\\overset{}{\\rightarrow}R^{D}$。\n", "MindSpore在自动微分的过程中，会计算向量-雅可比积。为了计算映射$R^{D}\\overset{}{\\rightarrow}R^{D}$的完整雅可比矩阵，我们每次使用不同的单位向量逐行来计算它。" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(16, 16)\n", "[-3.2045446e-05 -1.3530695e-05 1.8671712e-05 -9.6547810e-05\n", " 5.9755850e-05 -5.1343523e-05 1.3528993e-05 -4.6988782e-05\n", " -4.5517798e-05 -6.1188715e-05 -1.6264191e-04 5.5033437e-05\n", " -4.3497541e-05 2.2357668e-05 -1.3188722e-04 -3.0677278e-05]\n" ] } ], "source": [ "def partial_forecast(x):\n", " return ops.dense(x, weight, bias).tanh()\n", "\n", "_, vjp_fn = vjp(partial_forecast, x)\n", "\n", "def compute_jac_matrix(unit_vectors):\n", " jacobian_rows = [vjp_fn(vec)[0] for vec in unit_vectors]\n", " return ops.stack(jacobian_rows)\n", "\n", "\n", "unit_vectors = ops.eye(D)\n", "jacobian = compute_jac_matrix(unit_vectors)\n", "print(jacobian.shape)\n", "print(jacobian[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "在`compute_jac_matrix`中，使用for循环逐行计算的方式计算雅可比矩阵，计算效率并不高。MindSpore提供`jacrev`来计算雅可比矩阵，`jacrev`的实现利用了`vmap`，`vmap`可以消除`compute_jac_matrix`中的for循环并向量化整个计算过程。`jacrev`的参数`grad_position`指定计算输出相对于哪个参数的雅可比矩阵。\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "from mindspore import jacrev\n", "jacrev_jacobian = jacrev(forecast, grad_position=2)(weight, bias, x)\n", "assert np.allclose(jacrev_jacobian.asnumpy(), jacobian.asnumpy())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "接下来对`compute_jac_matrix`和`jacrev`的性能进行对比。通常情况下，`jacrev`的性能更好，因为使用`vmap`进行向量化运算，会更充分地利用硬件，同时计算多组数据，降低计算开销，以获得更好的性能。\n", "\n", "让我们编写一个函数，在微秒量级上评估两种方法的性能。" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "compute_jac_matrix run 500 times, cost time 12942823.04868102 microseconds.\n", "jacrev run 500 times, cost time 909309.7001314163 microseconds.\n" ] } ], "source": [ "def perf_compution(func, run_times, *args, **kwargs):\n", " start_time = time.perf_counter()\n", " for _ in range(run_times):\n", " func(*args, **kwargs)\n", " end_time = time.perf_counter()\n", " cost_time = (end_time - start_time) * 1000000\n", " return cost_time\n", "\n", "\n", "run_times = 500\n", "xp = x.copy()\n", "compute_jac_matrix_cost_time = perf_compution(compute_jac_matrix, run_times, xp)\n", "jac_fn = jacrev(forecast, grad_position=2)\n", "jacrev_cost_time = perf_compution(jac_fn, run_times, weight, bias, x)\n", "print(f\"compute_jac_matrix run {run_times} times, cost time {compute_jac_matrix_cost_time} microseconds.\")\n", "print(f\"jacrev run {run_times} times, cost time {jacrev_cost_time} microseconds.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "分别运行`compute_jac_matrix`和`jacrev`500次，统计它们消耗的时间。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "下面计算，相较于使用`compute_jac_matrix`，使用`jacrev`计算雅可比矩阵，性能提升的百分比。" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Performance delta: 92.9744 percent improvement with jacrev. \n" ] } ], "source": [ "def perf_cmp(first, first_descriptor, second, second_descriptor):\n", " faster = second\n", " slower = first\n", " gain = (slower - faster) / slower\n", " if gain < 0:\n", " gain *= -1\n", " final_gain = gain*100\n", " print(f\" Performance delta: {final_gain:.4f} percent improvement with {second_descriptor}. \")\n", "\n", "perf_cmp(compute_jac_matrix_cost_time, \"for loop\", jacrev_cost_time, \"jacrev\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "此外，也可以通过指定`jacrev`的参数`grad_position`来计算输出相对于模型参数weight和bias的雅可比矩阵。" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(16, 16, 16)\n", "(16, 16)\n" ] } ], "source": [ "jacrev_weight, jacrev_bias = jacrev(forecast, grad_position=(0, 1))(weight, bias, x)\n", "print(jacrev_weight.shape)\n", "print(jacrev_bias.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 反向模式计算雅可比矩阵 vs 前向模式计算雅可比矩阵\n", "\n", "MindSpore提供了两个API来计算雅可比矩阵：分别是`jacrev`和`jacfwd`。\n", "\n", "* `jacrev`：使用反向模式自动微分。\n", "\n", "* `jacfwd`：使用前向模式自动微分。\n", "\n", "`jacfwd`和`jacrev`可以相互替换，但是它们在不同的场景下，性能表现不同。\n", "\n", "一般来说，如果需要计算函数$R^{n}\\overset{}{\\rightarrow}R^{m}$的雅可比矩阵，当该函数的输出向量的规模大于输入向量的规模时（即，m > n），`jacfwd`在性能方面表现得更好，否则，`jacrev`在性能方面表现得更好。\n", "\n", "下面对这个结论做一个非严谨的论证，在前向模式自动微分（计算雅可比-向量积）的过程中，是逐列计算雅可比矩阵的，在反向模式自动微分（计算向量-雅可比积）的过程中，是逐行计算雅可比矩阵的。假设待计算的雅可比矩阵的规模是m行，n列，如果m > n，我们推荐使用逐列计算雅可比矩阵的`jacfwd`，反之，如果m < n，我们推荐使用逐行计算雅可比矩阵的`jacrev`。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 黑塞矩阵\n", "\n", "在介绍MindSpore提供的计算黑塞矩阵的方法之前，首先对黑塞矩阵进行介绍。\n", "\n", "黑塞矩阵可以由梯度操作$\\nabla$和广度梯度操作$\\partial$的复合得到，即\n", "\n", "$$\\nabla \\circ \\partial: F_{n}^{1} \\longrightarrow F_{n}^{n} \\longrightarrow F_{n \\times n}^{n}$$\n", "\n", "将该复合操作用于f，得到，\n", "\n", "$$f \\longmapsto \\nabla f \\longmapsto J_{\\nabla f}$$\n", "\n", "可以得到黑塞矩阵，\n", "\n", "$$H_{f} = \\begin{bmatrix}\n", " \\frac{\\partial (\\nabla _{1}f)}{\\partial x_{1}} &\\frac{\\partial (\\nabla _{1}f)}{\\partial x_{2}} &\\dots &\\frac{\\partial (\\nabla _{1}f)}{\\partial x_{n}} \\\\\n", " \\frac{\\partial (\\nabla _{2}f)}{\\partial x_{1}} &\\frac{\\partial (\\nabla _{2}f)}{\\partial x_{2}} &\\dots &\\frac{\\partial (\\nabla _{2}f)}{\\partial x_{n}} \\\\\n", " \\vdots &\\vdots &\\ddots &\\vdots \\\\\n", " \\frac{\\partial (\\nabla _{n}f)}{\\partial x_{1}} &\\frac{\\partial (\\nabla _{n}f)}{\\partial x_{2}} &\\dots &\\frac{\\partial (\\nabla _{n}f)}{\\partial x_{n}}\n", "\\end{bmatrix} = \\begin{bmatrix}\n", " \\frac{\\partial ^2 f}{\\partial x_{1}^{2}} &\\frac{\\partial ^2 f}{\\partial x_{2} \\partial x_{1}} &\\dots &\\frac{\\partial ^2 f}{\\partial x_{n} \\partial x_{1}} \\\\\n", " \\frac{\\partial ^2 f}{\\partial x_{1} \\partial x_{2}} &\\frac{\\partial ^2 f}{\\partial x_{2}^{2}} &\\dots &\\frac{\\partial ^2 f}{\\partial x_{n} \\partial x_{2}} \\\\\n", " \\vdots &\\vdots &\\ddots &\\vdots \\\\\n", " \\frac{\\partial ^2 f}{\\partial x_{1} \\partial x_{n}} &\\frac{\\partial ^2 f}{\\partial x_{2} \\partial x_{n}} &\\dots &\\frac{\\partial ^2 f}{\\partial x_{n}^{2}}\n", "\\end{bmatrix}$$\n", "\n", "易见，黑塞矩阵是一个实对称矩阵。\n", "\n", "黑塞矩阵的应用：利用黑塞矩阵，我们可以探索神经网络在某点处的曲率，为训练是否收敛提供数值依据。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 计算黑塞矩阵\n", "\n", "在MindSpore中，我们可以通过`jacfwd`和`jacrev`的任意组合来计算黑塞矩阵。" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "Din = 32\n", "Dout = 16\n", "weight = ops.randn(Dout, Din)\n", "bias = ops.randn(Dout)\n", "x = ops.randn(Din)\n", "\n", "hess1 = jacfwd(jacfwd(forecast, grad_position=2), grad_position=2)(weight, bias, x)\n", "hess2 = jacfwd(jacrev(forecast, grad_position=2), grad_position=2)(weight, bias, x)\n", "hess3 = jacrev(jacfwd(forecast, grad_position=2), grad_position=2)(weight, bias, x)\n", "hess4 = jacrev(jacrev(forecast, grad_position=2), grad_position=2)(weight, bias, x)\n", "\n", "np.allclose(hess1.asnumpy(), hess2.asnumpy())\n", "np.allclose(hess2.asnumpy(), hess3.asnumpy())\n", "np.allclose(hess3.asnumpy(), hess4.asnumpy())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 计算批量雅可比矩阵和批量黑塞矩阵\n", "\n", "在上面给出的示例中，我们都是计算单一的输出向量关于单一的输入向量的雅可比矩阵。在某些情况下，你可能想计算一批量的输出向量关于一批量的输入向量的雅可比矩阵，或者换句话说，给出一批量的输入向量，其shape为(b, n)，函数的映射关系是$R^{n}\\overset{}{\\rightarrow}R^{m}$，我们期望得到一批量的雅可比矩阵，其shape为(b, m, n)。\n", "\n", "我们可以使用`vmap`计算批量雅可比矩阵。" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(64, 33, 31)\n" ] } ], "source": [ "batch_size = 64\n", "Din = 31\n", "Dout = 33\n", "\n", "weight = ops.randn(Dout, Din)\n", "bias = ops.randn(Dout)\n", "x = ops.randn(batch_size, Din)\n", "\n", "compute_batch_jacobian = vmap(jacrev(forecast, grad_position=2), in_axes=(None, None, 0))\n", "batch_jacobian = compute_batch_jacobian(weight, bias, x)\n", "print(batch_jacobian.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "计算批量黑塞矩阵的方式与计算批量雅可比矩阵的方法类似，也可以使用`vmap`来计算批量黑塞矩阵，" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(64, 33, 31, 31)\n" ] } ], "source": [ "hessian = jacrev(jacrev(forecast, grad_position=2), grad_position=2)\n", "compute_batch_hessian = vmap(hessian, in_axes=(None, None, 0))\n", "batch_hessian = compute_batch_hessian(weight, bias, x)\n", "print(batch_hessian.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 计算黑塞-向量积\n", "\n", "计算黑塞-向量积（Hessian-vector product, hvp）的最直接的方法计算一个完整的黑塞矩阵，并将其与向量进行点积运算。但MindSpore提供了更好的方法，使得不需要计算一个完整的黑塞矩阵，便可以计算黑塞-向量积。下面我们介绍计算黑塞-向量积的两种方法。\n", "\n", "* 将反向模式自动微分与反向模式自动微分组合。\n", "\n", "* 将反向模式自动微分与前向模式自动微分组合。\n", "\n", "下面先介绍，在MindSpore中，如何使用反向模式自动微分与前向模式自动微分组合的方式计算黑塞-向量积，" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(128,)\n" ] } ], "source": [ "def hvp_revfwd(f, inputs, vector):\n", " return jvp(grad(f), inputs, vector)[1]\n", "\n", "def f(x):\n", " return x.sin().sum()\n", "\n", "inputs = ops.randn(128)\n", "vector = ops.randn(128)\n", "\n", "result_hvp_revfwd = hvp_revfwd(f, inputs, vector)\n", "print(result_hvp_revfwd.shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "如果前向自动微分不能满足要求，我们可以使用反向模式自动微分与反向模式自动微分组合的方式来计算黑塞-向量积，" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(128,)\n" ] } ], "source": [ "def hvp_revrev(f, inputs, vector):\n", " _, vjp_fn = vjp(grad(f), *inputs)\n", " return vjp_fn(*vector)\n", "\n", "result_hvp_revrev = hvp_revrev(f, (inputs,), (vector,))\n", "print(result_hvp_revrev[0].shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "使用上面两种方法计算黑塞-向量积得到的结果是一样的。" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "assert np.allclose(result_hvp_revfwd.asnumpy(), result_hvp_revrev[0].asnumpy())" ] } ], "metadata": { "kernelspec": { "display_name": "MindSpore", "language": "python", "name": "mindspore" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.10" } }, "nbformat": 4, "nbformat_minor": 5 }