diff --git a/.translate/state/newton_method.md.yml b/.translate/state/newton_method.md.yml new file mode 100644 index 00000000..26d1bc6d --- /dev/null +++ b/.translate/state/newton_method.md.yml @@ -0,0 +1,6 @@ +source-sha: 6c93e967f49a9104f11aab2447bd4ee4f0995b1c +synced-at: "2026-07-18" +model: claude-sonnet-5 +mode: RESYNC +section-count: 5 +tool-version: 0.17.0 diff --git a/lectures/newton_method.md b/lectures/newton_method.md index 7e57d728..84619c49 100644 --- a/lectures/newton_method.md +++ b/lectures/newton_method.md @@ -9,6 +9,26 @@ kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 +translation: + title: 使用牛顿法求解经济模型 + headings: + Overview: 概述 + Fixed point computation using Newton's method: 用牛顿法计算不动点 + Fixed point computation using Newton's method::The Solow model: 索洛模型 + Fixed point computation using Newton's method::Implementation: 实现 + Fixed point computation using Newton's method::Implementation::Successive approximation: 连续近似法 + Fixed point computation using Newton's method::Implementation::Newton's method: 牛顿法 + Root-Finding in one dimension: 一维求根 + Root-Finding in one dimension::Newton's method for zeros: 牛顿法求零点 + Root-Finding in one dimension::Application to finding fixed points: 在寻找不动点中的应用 + Multivariate Newton's method: 多元牛顿法 + Multivariate Newton's method::A two-goods market equilibrium: 双商品市场均衡 + Multivariate Newton's method::A two-goods market equilibrium::A graphical exploration: 图形化探索 + Multivariate Newton's method::A two-goods market equilibrium::Using a multidimensional root finder: 使用多维根查找器 + Multivariate Newton's method::A two-goods market equilibrium::Adding gradient information: 添加梯度信息 + Multivariate Newton's method::A two-goods market equilibrium::Using Newton's method: 使用牛顿法 + Multivariate Newton's method::A high-dimensional problem: 高维问题 + Exercises: 练习 --- (newton_method)= @@ -24,12 +44,11 @@ kernelspec: # 使用牛顿法求解经济模型 -```{contents} 目录 -:depth: 2 +```{include} _admonition/gpu.md ``` -```{seealso} -**GPU加速:** 本讲座有一个使用[jax](https://jax.readthedocs.io)在GPU上运行的版本，[点击这里查看](https://jax.quantecon.org/newtons_method.html) +```{contents} 目录 +:depth: 2 ``` ## 概述 @@ -48,23 +67,15 @@ kernelspec: 本讲将在一维和多维环境中应用牛顿法来解决不动点和零点计算问题。 -牛顿法的基本思路是: +* 在寻找函数 $f$ 的不动点时，牛顿法通过求解对函数 $f$ 的线性近似的不动点，来更新已有的不动点猜测值。 -* 对于不动点问题，通过对函数 $f$ 进行线性近似来寻找不动点。 - -* 对于零点问题，通过求解函数 $f$ 的线性近似的零点，不断更新当前的估计值，直到收敛到真实的零点。 +* 在寻找函数 $f$ 的零点时，牛顿法通过求解对函数 $f$ 的线性近似的零点，来更新已有的猜测值。 为了建立直观认识，我们首先考虑一个简单的一维不动点问题，其中我们已知解，并使用连续近似和牛顿法来求解。 然后我们将牛顿法应用到多维环境中，求解多种商品的市场均衡问题。 -最后，我们将使用 [`autograd`](https://github.com/HIPS/autograd) 包提供的自动微分功能来处理一个高维均衡问题。 - -```{code-cell} ipython3 -:tags: [hide-output] - -!pip install autograd -``` +在本讲的最后，我们将利用 [`jax`](https://docs.jax.dev/en/latest/_autosummary/jax.grad.html) 中自动微分的能力来求解一个非常高维的均衡问题。 我们在本讲中使用以下导入语句 @@ -74,14 +85,13 @@ import matplotlib as mpl FONTPATH = "fonts/SourceHanSerifSC-SemiBold.otf" mpl.font_manager.fontManager.addfont(FONTPATH) plt.rcParams['font.family'] = ['Source Han Serif SC'] - -from collections import namedtuple +from typing import NamedTuple from scipy.optimize import root -from autograd import jacobian -# 经过简单封装的numpy，以支持自动微分 -import autograd.numpy as np +import jax.numpy as jnp +import jax -plt.rcParams["figure.figsize"] = (10, 5.7) +# 启用64位精度 +jax.config.update("jax_enable_x64", True) ``` ## 用牛顿法计算不动点 @@ -117,21 +127,27 @@ plt.rcParams["figure.figsize"] = (10, 5.7) 用纸笔解方程 $g(k)=k$，你可以验证 -$$ k^* = \left(\frac{s A}{δ}\right)^{1/(1 - α)} $$ +$$ +k^* = \left(\frac{s A}{\delta}\right)^{1/(1 - \alpha)} +$$ ### 实现 -让我们使用 [`namedtuple`](https://docs.python.org/3/library/collections.html#collections.namedtuple) 来存储我们的参数，这有助于保持代码的整洁和简洁。 +让我们使用 [`NamedTuple`](https://typing.python.org/en/latest/spec/namedtuples.html) 来存储我们的参数，这有助于保持代码的整洁和简洁。 ```{code-cell} ipython3 -SolowParameters = namedtuple("SolowParameters", ('A', 's', 'α', 'δ')) +class SolowParameters(NamedTuple): + A: float + s: float + α: float + δ: float ``` -此函数创建一个带有默认参数值的适当的`namedtuple`。 +此函数创建一个带有默认参数值的适当的 `SolowParameters`。 ```{code-cell} ipython3 def create_solow_params(A=2.0, s=0.3, α=0.3, δ=0.4): - "创建一个带有默认值的索洛模型参数化。" + """Creates a Solow model parameterization with default values.""" return SolowParameters(A=A, s=s, α=α, δ=δ) ``` @@ -141,35 +157,38 @@ def create_solow_params(A=2.0, s=0.3, α=0.3, δ=0.4): def g(k, params): A, s, α, δ = params return A * s * k**α + (1 - δ) * k - + + def exact_fixed_point(params): A, s, α, δ = params - return ((s * A) / δ)**(1/(1 - α)) + return ((s * A) / δ) ** (1 / (1 - α)) ``` 这是一个用于绘制45度动态图的函数。 ```{code-cell} ipython3 def plot_45(params, ax, fontsize=14): - + k_min, k_max = 0.0, 3.0 - k_grid = np.linspace(k_min, k_max, 1200) + k_grid = jnp.linspace(k_min, k_max, 1200) # 绘制函数 lb = r"$g(k) = sAk^{\alpha} + (1 - \delta)k$" - ax.plot(k_grid, g(k_grid, params), lw=2, alpha=0.6, label=lb) + ax.plot(k_grid, g(k_grid, params), lw=2, alpha=0.6, label=lb) ax.plot(k_grid, k_grid, "k--", lw=1, alpha=0.7, label="45") # 显示并标注固定点 kstar = exact_fixed_point(params) fps = (kstar,) ax.plot(fps, fps, "go", ms=10, alpha=0.6) - ax.annotate(r"$k^* = (sA / \delta)^{\frac{1}{1-\alpha}}$", - xy=(kstar, kstar), - xycoords="data", - xytext=(20, -20), - textcoords="offset points", - fontsize=fontsize) + ax.annotate( + r"$k^* = (sA / \delta)^{\frac{1}{1-\alpha}}$", + xy=(kstar, kstar), + xycoords="data", + xytext=(20, -20), + textcoords="offset points", + fontsize=fontsize, + ) ax.legend(loc="upper left", frameon=False, fontsize=fontsize) @@ -209,7 +228,7 @@ plt.show() ```{code-cell} ipython3 def compute_iterates(k_0, f, params, n=25): - "计算由任意函数f生成的长度为n的时间序列。" + """Compute time series of length n generated by function f.""" k = k_0 k_iterates = [] for t in range(n): @@ -225,8 +244,8 @@ k_series = compute_iterates(k_0, g, params) k_star = exact_fixed_point(params) fig, ax = plt.subplots() -ax.plot(k_series, 'o') -ax.plot([k_star] * len(k_series), 'k--') +ax.plot(k_series, "o") +ax.plot([k_star] * len(k_series), "k--") ax.set_ylim(0, 3) plt.show() ``` @@ -241,6 +260,8 @@ k_star_approx 这接近真实值。 +(solved_k)= + ```{code-cell} ipython3 k_star ``` @@ -254,13 +275,13 @@ k_star ```{math} :label: motivation -\hat g(x) \approx g(x_0)+g'(x_0)(x-x_0) +\hat g(x) \approx g(x_0) + g'(x_0)(x - x_0) ``` 我们通过计算满足以下等式的$x_1$来求解$\hat g$的不动点： $$ -x_1=\frac{g(x_0)-g'(x_0) x_0}{1-g'(x_0)} +x_1 = \frac{g(x_0) - g'(x_0) x_0}{1 - g'(x_0)} $$ 推广上述过程，牛顿不动点法的迭代公式为： @@ -277,7 +298,7 @@ x_{t+1} = \frac{g(x_t) - g'(x_t) x_t}{ 1 - g'(x_t) }, ```{math} :label: newton_method2 -g'(k) = \alpha s A k^{\alpha-1} + (1-\delta) +g'(k) = \alpha s A k^{\alpha - 1} + (1 - \delta) ``` @@ -286,7 +307,7 @@ g'(k) = \alpha s A k^{\alpha-1} + (1-\delta) ```{code-cell} ipython3 def Dg(k, params): A, s, α, δ = params - return α * A * s * k**(α-1) + (1 - δ) + return α * A * s * k ** (α - 1) + (1 - δ) ``` 下面的函数 $q$ 表示 [](newtons_method)。 @@ -299,11 +320,13 @@ def q(k, params): 现在让我们绘制一些轨迹。 ```{code-cell} ipython3 -def plot_trajectories(params, - k0_a=0.8, # 第一个初始条件 - k0_b=3.1, # 第二个初始条件 - n=20, # 时间序列长度 - fs=14): # 字体大小 +def plot_trajectories( + params, + k0_a=0.8, # 第一个初始条件 + k0_b=3.1, # 第二个初始条件 + n=20, # 时间序列长度 + fs=14, # 字体大小 +): fig, axes = plt.subplots(2, 1, figsize=(10, 6)) ax1, ax2 = axes @@ -321,13 +344,13 @@ def plot_trajectories(params, ax2.plot(ks4, "-o", label="牛顿法") for ax in axes: - ax.plot(k_star * np.ones(n), "k--") + ax.plot(k_star * jnp.ones(n), "k--") ax.legend(fontsize=fs, frameon=False) ax.set_ylim(0.6, 3.2) ax.set_yticks((k_star,)) ax.set_yticklabels(("$k^*$",), fontsize=fs) - ax.set_xticks(np.linspace(0, 19, 20)) - + ax.set_xticks(jnp.linspace(0, 19, 20)) + plt.show() ``` @@ -379,11 +402,15 @@ x_{t+1} = x_t - \frac{ f(x_t) }{ f'(x_t) }, 以下代码实现了迭代公式 [](oneD-newton) +(first_newton_attempt)= + ```{code-cell} ipython3 -def newton(f, Df, x_0, tol=1e-7, max_iter=100_000): +def newton(f, x_0, tol=1e-7, max_iter=100_000): x = x_0 + Df = jax.grad(f) # 实现零点查找公式 + @jax.jit def q(x): return x - f(x) / Df(x) @@ -391,13 +418,13 @@ def newton(f, Df, x_0, tol=1e-7, max_iter=100_000): n = 0 while error > tol: n += 1 - if(n > max_iter): + if n > max_iter: raise Exception('达到最大迭代次数但未收敛') y = q(x) - error = np.abs(x - y) + error = jnp.abs(x - y) x = y print(f'迭代 {n}, 误差 = {error:.5f}') - return x + return x.item() ``` 许多库都实现了一维牛顿法，包括SciPy，所以这里的代码仅作说明用途。 @@ -416,9 +443,7 @@ def newton(f, Df, x_0, tol=1e-7, max_iter=100_000): ```{code-cell} ipython3 params = create_solow_params() -k_star_approx_newton = newton(f=lambda x: g(x, params) - x, - Df=lambda x: Dg(x, params) - 1, - x_0=0.8) +k_star_approx_newton = newton(f = lambda x: g(x, params) - x, x_0=0.8) ``` ```{code-cell} ipython3 @@ -437,6 +462,7 @@ k_star_approx_newton 我们将看到使用牛顿法时能获得显著的性能提升。 +(two_goods_market)= ### 双商品市场均衡 让我们从计算双商品问题的市场均衡开始。 @@ -470,20 +496,20 @@ $$ 我们设定 $$ -A = \begin{pmatrix} +A = \begin{bmatrix} a_{00} & a_{01} \\ a_{10} & a_{11} - \end{pmatrix}, + \end{bmatrix}, \qquad - b = \begin{pmatrix} + b = \begin{bmatrix} b_0 \\ b_1 - \end{pmatrix} + \end{bmatrix} \qquad \text{和} \qquad - c = \begin{pmatrix} + c = \begin{bmatrix} c_0 \\ c_1 - \end{pmatrix} + \end{bmatrix} $$ 用于这个特定问题。 @@ -496,55 +522,66 @@ $$ $$ e(p) = - \begin{pmatrix} + \begin{bmatrix} e_0(p) \\ e_1(p) - \end{pmatrix} + \end{bmatrix} $$ 下面的函数计算给定参数的超额需求 ```{code-cell} ipython3 +@jax.jit def e(p, A, b, c): - return np.exp(- A @ p) + c - b * np.sqrt(p) + return jnp.exp(-A @ p) + c - b * jnp.sqrt(p) ``` 我们的默认参数值将是 $$ -A = \begin{pmatrix} +A = \begin{bmatrix} 0.5 & 0.4 \\ 0.8 & 0.2 - \end{pmatrix}, + \end{bmatrix}, \qquad - b = \begin{pmatrix} + b = \begin{bmatrix} 1 \\ 1 - \end{pmatrix} + \end{bmatrix} \qquad \text{和} \qquad - c = \begin{pmatrix} + c = \begin{bmatrix} 1 \\ 1 - \end{pmatrix} + \end{bmatrix} $$ ```{code-cell} ipython3 -A = np.array([ - [0.5, 0.4], - [0.8, 0.2] -]) -b = np.ones(2) -c = np.ones(2) +A = jnp.array([[0.5, 0.4], [0.8, 0.2]]) +b = jnp.ones(2) +c = jnp.ones(2) ``` 在价格水平 $p = (1, 0.5)$ 时，超额需求为 ```{code-cell} ipython3 -ex_demand = e((1.0, 0.5), A, b, c) +p = jnp.array([1, 0.5]) +ex_demand = e(p, A, b, c) -print(f'商品0的超额需求为 {ex_demand[0]:.3f} \n' - f'商品1的超额需求为 {ex_demand[1]:.3f}') +print( + f"商品0的超额需求为 {ex_demand[0]:.3f} \n" + f"商品1的超额需求为 {ex_demand[1]:.3f}" +) +``` + +为了提高计算效率，我们将使用 [`jax.vmap`](https://docs.jax.dev/en/latest/_autosummary/jax.vmap.html) 提供的向量化能力。这比使用python循环快得多。 + +```{code-cell} ipython3 +# 在p的第一个轴上创建向量化。 +e_vectorized_p_1 = jax.vmap(e, in_axes=(0, None, None, None)) + +# 在p的第二个轴上创建向量化。 +e_vectorized = jax.vmap(e_vectorized_p_1, in_axes=(0, None, None, None)) ``` 接下来我们在$(p_0, p_1)$值的网格上绘制两个函数$e_0$和$e_1$的等高线图和曲面图。 @@ -553,14 +590,17 @@ print(f'商品0的超额需求为 {ex_demand[0]:.3f} \n' ```{code-cell} ipython3 def plot_excess_demand(ax, good=0, grid_size=100, grid_max=4, surface=True): + p_grid = jnp.linspace(0, grid_max, grid_size) - # 创建一个100x100的网格 - p_grid = np.linspace(0, grid_max, grid_size) - z = np.empty((100, 100)) + # 为p_1和p_2的所有组合创建网格 + P1, P2 = jnp.meshgrid(p_grid, p_grid, indexing="ij") - for i, p_1 in enumerate(p_grid): - for j, p_2 in enumerate(p_grid): - z[i, j] = e((p_1, p_2), A, b, c)[good] + # 堆叠成形状为(grid_size, grid_size, 2)的数组 + P = jnp.stack([P1, P2], axis=-1) + + # 使用向量化函数一次性计算所有值 + z_full = e_vectorized(P, A, b, c) + z = z_full[:, :, good] if surface: cs1 = ax.contourf(p_grid, p_grid, z.T, alpha=0.5) @@ -611,14 +651,14 @@ plt.show() 我们以 $p = (1, 1)$ 作为初始猜测值。 ```{code-cell} ipython3 -init_p = np.ones(2) +init_p = jnp.ones(2) ``` 这个算法使用[改进的Powell方法](https://docs.scipy.org/doc/scipy/reference/optimize.root-hybr.html#optimize-root-hybr)来寻找零点 ```{code-cell} ipython3 %%time -solution = root(lambda p: e(p, A, b, c), init_p, method='hybr') +solution = root(lambda p: e(p, A, b, c), init_p, method="hybr") ``` 这是得到的值 @@ -631,7 +671,8 @@ p 这个结果看起来和我们从图中观察到的猜测很接近。我们可以把它代回到 $e$ 中验证 $e(p) \approx 0$ ```{code-cell} ipython3 -np.max(np.abs(e(p, A, b, c))) +e_p = jnp.max(jnp.abs(e(p, A, b, c))) +e_p.item() ``` 这确实是一个很小的误差。 @@ -645,10 +686,10 @@ np.max(np.abs(e(p, A, b, c))) $$ J(p) = - \begin{pmatrix} + \begin{bmatrix} \frac{\partial e_0}{\partial p_0}(p) & \frac{\partial e_0}{\partial p_1}(p) \\ \frac{\partial e_1}{\partial p_0}(p) & \frac{\partial e_1}{\partial p_1}(p) - \end{pmatrix} + \end{bmatrix} $$ ```{code-cell} ipython3 @@ -656,28 +697,30 @@ def jacobian_e(p, A, b, c): p_0, p_1 = p a_00, a_01 = A[0, :] a_10, a_11 = A[1, :] - j_00 = -a_00 * np.exp(-a_00 * p_0) - (b[0]/2) * p_0**(-1/2) - j_01 = -a_01 * np.exp(-a_01 * p_1) - j_10 = -a_10 * np.exp(-a_10 * p_0) - j_11 = -a_11 * np.exp(-a_11 * p_1) - (b[1]/2) * p_1**(-1/2) - J = [[j_00, j_01], - [j_10, j_11]] - return np.array(J) + j_00 = -a_00 * jnp.exp(-a_00 * p_0) - (b[0] / 2) * p_0 ** (-1 / 2) + j_01 = -a_01 * jnp.exp(-a_01 * p_1) + j_10 = -a_10 * jnp.exp(-a_10 * p_0) + j_11 = -a_11 * jnp.exp(-a_11 * p_1) - (b[1] / 2) * p_1 ** (-1 / 2) + J = [[j_00, j_01], [j_10, j_11]] + return jnp.array(J) ``` ```{code-cell} ipython3 %%time -solution = root(lambda p: e(p, A, b, c), - init_p, - jac=lambda p: jacobian_e(p, A, b, c), - method='hybr') +solution = root( + lambda p: e(p, A, b, c), + init_p, + jac = lambda p: jacobian_e(p, A, b, c), + method="hybr", +) ``` 现在的解更加精确了（尽管在这个低维问题中，差异非常小）： ```{code-cell} ipython3 p = solution.x -np.max(np.abs(e(p, A, b, c))) +e_p = jnp.max(jnp.abs(e(p, A, b, c))) +e_p.item() ``` #### 使用牛顿法 @@ -696,35 +739,35 @@ p_{n+1} = p_n - J_e(p_n)^{-1} e(p_n) 迭代从价格向量$p_0$的某个初始猜测开始。 -在这里，我们不手动编写雅可比矩阵，而是使用`autograd`库中的`jacobian()`函数来自动求导并计算雅可比矩阵。 +在这里，我们不手动编写雅可比矩阵，而是使用`jax`库中的`jacobian()`函数来自动求导并计算雅可比矩阵。 只需稍作修改，我们就可以将[我们之前的尝试](first_newton_attempt)推广到多维问题 ```{code-cell} ipython3 def newton(f, x_0, tol=1e-5, max_iter=10): x = x_0 - q = lambda x: x - np.linalg.solve(jacobian(f)(x), f(x)) + f_jac = jax.jacobian(f) + + @jax.jit + def q(x): + return x - jnp.linalg.solve(f_jac(x), f(x)) + error = tol + 1 n = 0 while error > tol: - n+=1 - if(n > max_iter): - raise Exception('Max iteration reached without convergence') + n += 1 + if n > max_iter: + raise Exception("Max iteration reached without convergence") y = q(x) - if(any(np.isnan(y))): - raise Exception('Solution not found with NaN generated') - error = np.linalg.norm(x - y) + if any(jnp.isnan(y)): + raise Exception("Solution not found with NaN generated") + error = jnp.linalg.norm(x - y) x = y - print(f'iteration {n}, error = {error:.5f}') - print('\n' + f'Result = {x} \n') + print(f'迭代 {n}, 误差 = {error:.5f}') + print("\n" + f"Result = {x} \n") return x ``` -```{code-cell} ipython3 -def e(p, A, b, c): - return np.exp(- np.dot(A, p)) + c - b * np.sqrt(p) -``` - 我们发现算法在4步内终止 ```{code-cell} ipython3 @@ -733,38 +776,41 @@ p = newton(lambda p: e(p, A, b, c), init_p) ``` ```{code-cell} ipython3 -np.max(np.abs(e(p, A, b, c))) +e_p = jnp.max(jnp.abs(e(p, A, b, c))) +e_p.item() ``` 结果非常准确。 +在开销较大的情况下，速度并不比经过优化的`scipy`函数更快。 + + ### 高维问题 我们的下一步是研究一个有3,000种商品的大型市场。 -使用GPU加速线性代数和自动微分的JAX版本可在[此处](https://jax.quantecon.org/newtons_method.html#application)获取 - 超额需求函数基本相同，但现在矩阵 $A$ 是 $3000 \times 3000$ 的，参数向量 $b$ 和 $c$ 是 $3000 \times 1$ 的。 ```{code-cell} ipython3 dim = 3000 -np.random.seed(123) -# 创建随机矩阵A并将行归一化使其和为1 -A = np.random.rand(dim, dim) -A = np.asarray(A) -s = np.sum(A, axis=0) +# 创建JAX随机密钥 +key = jax.random.PRNGKey(0) + +# 创建随机矩阵A并将列归一化使其和为1 +A = jax.random.uniform(key, (dim, dim)) +s = jnp.sum(A, axis=0) A = A / s # 设置b和c -b = np.ones(dim) -c = np.ones(dim) +b = jnp.ones(dim) +c = jnp.ones(dim) ``` 这是我们的初始条件 ```{code-cell} ipython3 -init_p = np.ones(dim) +init_p = jnp.ones(dim) ``` ```{code-cell} ipython3 @@ -773,23 +819,27 @@ p = newton(lambda p: e(p, A, b, c), init_p) ``` ```{code-cell} ipython3 -np.max(np.abs(e(p, A, b, c))) +e_p = jnp.max(jnp.abs(e(p, A, b, c))) +e_p.item() ``` 在相同的容差条件下，我们比较牛顿法与SciPy的`root`函数的运行时间和精确度 ```{code-cell} ipython3 %%time -solution = root(lambda p: e(p, A, b, c), - init_p, - jac=lambda p: jacobian(e)(p, A, b, c), - method='hybr', - tol=1e-5) +solution = root( + lambda p: e(p, A, b, c), + init_p, + jac = lambda p: jax.jacobian(e)(p, A, b, c), + method="hybr", + tol=1e-5, +) ``` ```{code-cell} ipython3 p = solution.x -np.max(np.abs(e(p, A, b, c))) +e_p = jnp.max(jnp.abs(e(p, A, b, c))) +e_p.item() ``` ## 练习 @@ -801,13 +851,13 @@ np.max(np.abs(e(p, A, b, c))) 考虑索洛固定点问题的三维扩展，其中 $$ -A = \begin{pmatrix} +A = \begin{bmatrix} 2 & 3 & 3 \\ 2 & 4 & 2 \\ 1 & 5 & 1 \\ - \end{pmatrix}, + \end{bmatrix}, \quad -s = 0.2, \quad α = 0.5, \quad δ = 0.8 +s = 0.2, \quad \alpha = 0.5, \quad \delta = 0.8 $$ 和之前一样，运动方程为 @@ -829,23 +879,22 @@ $$ \end{aligned} $$ -````{hint} +````{hint} :class: dropdown -- 固定点的计算等价于计算满足 $f(k^*) - k^* = 0$ 的 $k^*$。 +- 固定点的计算等价于计算满足 $g(k^*) - k^* = 0$ 的 $k^*$。 - 如果你对你的解决方案不确定，可以从已解决的示例开始： ```{math} -A = \begin{pmatrix} +A = \begin{bmatrix} 2 & 0 & 0 \\ 0 & 2 & 0 \\ - -0 & 0 & 2 \\ - \end{pmatrix} + 0 & 0 & 2 \\ + \end{bmatrix} ``` -其中 $s = 0.3$、$α = 0.3$ 和 $δ = 0.4$，初始值为： +其中 $s = 0.3$、$\alpha = 0.3$ 和 $\delta = 0.4$，初始值为： ```{math} k_0 = (1, 1, 1) @@ -865,24 +914,21 @@ k_0 = (1, 1, 1) 让我们首先定义这个问题的参数 ```{code-cell} ipython3 -A = np.array([[2.0, 3.0, 3.0], - [2.0, 4.0, 2.0], - [1.0, 5.0, 1.0]]) +A = jnp.array([[2.0, 3.0, 3.0], [2.0, 4.0, 2.0], [1.0, 5.0, 1.0]]) s = 0.2 α = 0.5 δ = 0.8 -initLs = [np.ones(3), - np.array([3.0, 5.0, 5.0]), - np.repeat(50.0, 3)] +initLs = [jnp.ones(3), jnp.array([3.0, 5.0, 5.0]), jnp.repeat(50.0, 3)] ``` 然后定义[资本运动定律](motion_law)的多元版本 ```{code-cell} ipython3 +@jax.jit def multivariate_solow(k, A=A, s=s, α=α, δ=δ): - return (s * np.dot(A, k**α) + (1 - δ) * k) + return s * jnp.dot(A, k**α) + (1 - δ) * k ``` 让我们遍历每个初始值并查看输出结果 @@ -890,7 +936,7 @@ def multivariate_solow(k, A=A, s=s, α=α, δ=δ): ```{code-cell} ipython3 attempt = 1 for init in initLs: - print(f'尝试 {attempt}: 初始值为 {init} \n') + print(f'Attempt {attempt}: Starting value is {init} \n') %time k = newton(lambda k: multivariate_solow(k) - k, \ init) print('-'*64) @@ -912,27 +958,29 @@ multivariate_solow(k) - k 我们也可以在已知解上测试我们的结果 ```{code-cell} ipython3 -A = np.array([[2.0, 0.0, 0.0], - [0.0, 2.0, 0.0], - [0.0, 0.0, 2.0]]) +A = jnp.array([[2.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]]) s = 0.3 α = 0.3 δ = 0.4 -init = np.repeat(1.0, 3) +init = jnp.repeat(1.0, 3) +``` +```{code-cell} ipython3 +%%time -%time k = newton(lambda k: multivariate_solow(k, A=A, s=s, α=α, δ=δ) - k, \ - init) +k = newton(lambda k: multivariate_solow(k, A=A, s=s, α=α, δ=δ) - k, init) ``` 结果与真实值非常接近，但仍有细微差异。 ```{code-cell} ipython3 -%time k = newton(lambda k: multivariate_solow(k, A=A, s=s, α=α, δ=δ) - k, \ - init,\ - tol=1e-7) +%%time + +k = newton( + lambda k: multivariate_solow(k, A=A, s=s, α=α, δ=δ) - k, init, tol=1e-7 +) ``` 我们可以看到它正在朝着更精确的解迈进。 @@ -950,29 +998,28 @@ init = np.repeat(1.0, 3) 让我们定义一个具有以下默认值的三商品问题： $$ -A = \begin{pmatrix} +A = \begin{bmatrix} 0.2 & 0.1 & 0.7 \\ 0.3 & 0.2 & 0.5 \\ 0.1 & 0.8 & 0.1 \\ - \end{pmatrix}, + \end{bmatrix}, \qquad -b = \begin{pmatrix} +b = \begin{bmatrix} 1 \\ 1 \\ 1 - \end{pmatrix} + \end{bmatrix} \qquad \text{和} \qquad -c = \begin{pmatrix} +c = \begin{bmatrix} 1 \\ 1 \\ 1 - \end{pmatrix} + \end{bmatrix} $$ 对于这个练习，使用以下极端价格向量作为初始值： $$ - \begin{aligned} p1_{0} &= (5, 5, 5) \\ p2_{0} &= (1, 1, 1) \\ @@ -980,63 +1027,48 @@ $$ \end{aligned} $$ -将容差设置为$0.0$以获得更精确的输出。 +将容差设置为$1e-15$以获得更精确的输出。 ```{exercise-end} ``` ```{solution-start} newton_ex2 - :class: dropdown ``` 定义参数和初始值 ```{code-cell} ipython3 -A = np.array([ - [0.2, 0.1, 0.7], - [0.3, 0.2, 0.5], - [0.1, 0.8, 0.1] -]) - -b = np.array([1.0, 1.0, 1.0]) -c = np.array([1.0, 1.0, 1.0]) +A = jnp.array([[0.2, 0.1, 0.7], [0.3, 0.2, 0.5], [0.1, 0.8, 0.1]]) +b = jnp.array([1.0, 1.0, 1.0]) +c = jnp.array([1.0, 1.0, 1.0]) -initLs = [np.repeat(5.0, 3), - np.ones(3), - np.array([4.5, 0.1, 4.0])] +initLs = [jnp.repeat(5.0, 3), jnp.ones(3), jnp.array([4.5, 0.1, 4.0])] ``` 让我们检查每个初始猜测值并查看输出结果 ```{code-cell} ipython3 ---- -tags: [raises-exception] ---- +:tags: [raises-exception] -attempt = 1 -for init in initLs: - print(f'尝试 {attempt}: 初始值为 {init} \n') - %time p = newton(lambda p: e(p, A, b, c), \ - init, \ - tol=1e-15, \ - max_iter=15) - print('-'*64) - attempt += 1 +for attempt, init in enumerate(initLs, start=1): + print(f"Attempt {attempt}: Starting value is {init} \n") + %time p = newton(lambda p: e(p, A, b, c), init, tol=1e-15, max_iter=15) + print("-" * 64) ``` 我们可以发现牛顿法对某些初始值可能会失败。 有时可能需要尝试几个初始猜测值才能实现收敛。 -将结果代回公式中检验我们的结果 +将结果代回公式中检验我们的结果，使用收敛的第二个初始猜测值 ```{code-cell} ipython3 -e(p, A, b, c) +p_solution = newton(lambda p: e(p, A, b, c), initLs[1], tol=1e-15, max_iter=15) +e(p_solution, A, b, c) ``` 我们可以看到结果非常精确。 ```{solution-end} ``` -