优化：牛顿-拉夫森方法的几何解释
探索一种基本的数值优化技术，重点关注其几何解释
·发布于Towards Data Science ·阅读时间 8 分钟·2023 年 10 月 13 日
–
图片由Ansgar Scheffold拍摄，发布在Unsplash
梯度下降被广泛认为是基本的数值优化技术之一，而牛顿-拉夫森方法则在这一领域中具有显著地位。这种方法以其简单性、优雅性和计算能力而著称，值得深入探讨。
在本文中，我们的目标是阐明牛顿-拉夫森方法的几何原理。通过这种阐述，我们旨在为读者提供对其机制的直观理解，并消除与其数学基础相关的潜在复杂性。
随后，为了建立我们讨论的坚实数学框架，我们将深入探讨该方法的数学复杂性，并提供在 Python 编程语言中的实际实现。
在此演示之后，我们将区分牛顿-拉夫森方法的两个主要应用：根查找和优化。这一区分将明确方法在不同上下文中的实际应用。
最后，我们将对牛顿-拉夫森方法和梯度下降方法进行比较分析，提供对它们各自优缺点的见解。
如果你对数学概念感兴趣，并希望通过 Python 快速学习，请查看我的书籍：
[## 用 Python 揭示 70 个数学概念：通过 Python 探索数学的实用指南
购买《用 Python 揭示 70 个数学概念：通过 Python 探索数学的实用指南》请访问…
图形概述
查找根的迭代过程。图片由作者提供。
基本上，牛顿-拉夫森方法是一种迭代过程，旨在数值确定数学方程的根。其主要目标是确定一个值，记作 x_root，使方程 f(x_root) = 0 得以满足。需要注意的是，由于其数值性质，通过该方法获得的 x_root 值是一个近似值。
该方法的核心原理如下：从初始点 x_0 开始，目标是生成一系列值 x_i，这些值逐渐逼近所寻求的 x_root。为了实现这一点，在 (x_i, f(x_i)) 点附近通过线性化过程对函数进行局部近似。这种近似通过计算函数在点 (x_i, f(x_i)) 处的切线来完成。
对于这种局部近似，根的确定是直接的：它对应于该切线与 x 轴的交点。这个交点的横坐标 x_i+1 作为根 x_root 的更新近似值。
这个迭代过程会不断重复，直到 f(x) 的值足够接近零，表明已接近根。上面的图形作为该过程的示意图，展示了每次迭代时的切线及其与 x 轴的交点。
在这个具体的例子中，使用的基础函数是多项式 (x-6)(x-6)(x-6)，它在 x=6 处具有三重根。显然，经过几次迭代，大约 10 次，该方法会收敛到接近理论根 x_root=6 的点。
绘制该图形的代码如下：
显示牛顿-拉夫森方法的连续迭代。代码由作者提供。
按照我通常的做法，我在上述代码中依赖自动微分，特别是使用 JAX 库来有效计算 f 的导数。
方法的数学推导
在让你感受到牛顿-拉夫森方法的工作方式后，是时候看看这个方法如何从数学上推导出来了。
给定一个函数f(x)和一个起始点x_i，目的是在该点(x_i, f(x_i))对函数进行线性化。这正是我们在前一节中所做的，通过在该点绘制切线。
从形式上讲，这可以表示为f在x_i处的一阶泰勒展开：
一阶泰勒展开。作者提供的公式。
该公式给出了f在x_i附近的线性化公式，其中delta是变量。
给定此公式，按照上述几何方法计算delta，更新x_i，使得该线性化与 x 轴相交，即：
计算x_i+1。作者提供的公式。
牛顿-拉夫森方法用于根的识别
正如引言中提到的，牛顿-拉夫森方法可以用来解决两个主要问题：根的查找和函数优化。在本节中，我们将特别关注根的查找。
首先，让我们重新审视根的查找概念。在几何背景下，这本质上是识别一个函数f与 x 轴交点的横坐标x。这与本文初始图中的基本概念相符。从代数角度来看，这等同于确定x使得f(x) = 0。
因此，牛顿-拉夫森方法被设计用来解决这个具体问题：确定函数f的曲线与 x 轴交点的x值。
牛顿-拉夫森方法用于优化
如前所述，该方法主要解决的问题是根的查找。然而，在许多情况下，牛顿-拉夫森方法并不是用来寻找根，而是用于优化。换句话说，它用于识别x的值，使得函数f达到最大值或最小值点。
这怎么可能呢？最初，这可能看起来根的查找和优化是完全不同的挑战。然而，它们却是密切相关的。这种联系来源于一个事实，即函数在其斜率符号改变的点达到极值，这意味着其导数为零。因此，寻找函数的极值等同于识别其导数的根。
将牛顿-拉夫森方法应用于这种情况是相对简单的。该过程包括用其导数f’替换f，然后用其二阶导数f’’替换f’。
以下代码片段展示了如何调整该方法以最大化高度非线性函数-(-(x-6)**2 / 3.0 + log(6+x) + 2*sin(x))：
该代码通过应用牛顿-拉夫森方法来查找函数达到最大值的x，并绘制连续探索的点：
优化一个明显非线性的函数。作者提供的图像。
牛顿-拉夫森方法的优化几何解释
牛顿-拉夫森方法在优化中的适应可以从两种基本方式进行解释。第一种方法涉及将用于根查找的几何解释应用于函数 f 的导数 f'，可以在本文的早期部分找到详细信息。
第二种更具吸引力的观点认为，在优化的背景下，函数 f 不再由直线（即一阶多项式）局部近似，而是由二阶多项式近似。从几何角度来看，这意味着 f 的曲线由抛物线近似。在这种情况下，目标不是确定局部近似与 x 轴的交点（如同在根查找中），而是找到该函数的最小值。
正如俗话所说，“一图胜千言”，下图的视觉表示有助于有效地说明这个概念。
通过拟合抛物线进行优化。图由作者提供。
基本上，算法的工作原理如下：
-
选择一个点
-
在这个点拟合一个抛物线，即一个多项式
p(x)=ax**2+bx+c，使其通过同一个点，并且与f具有相同的一阶和二阶导数。 -
使用抛物线极值的横坐标
x作为新的起始点 -
迭代
从 Python 的角度来看，这给出：
牛顿-拉夫森方法的几何再解释。代码由作者提供。
请注意，这种牛顿-拉夫森方法的几何再解释给出的结果与之前的实现完全相同（在此情况下为7.4074383），并且收敛速度相同。只是稍微不够高效，因为我们需要计算多项式 p 的三个权重。
应用条件
牛顿-拉夫森方法是一个非常优雅且强大的方法。然而，在应用于优化时，它需要一些条件才能发挥作用。
首先，优化的函数必须是二次可微的。其次，需要从一个良好的近似开始，以避免局部最小值。
牛顿-拉夫森方法不同于梯度下降
另一个非常流行的优化方法是梯度下降。然而，它们在许多方面有所不同。
首先，牛顿-拉夫森方法在应用于优化时要求函数 f 是二次可微的。相对而言，梯度下降仅要求一阶导数。
其次，牛顿-拉夫森方法需要的调整较少，因为没有学习率需要配置。新点由抛物线的最小值自动给出。
第三，这种直接估计下一个点的方法，无需设置学习率，确保了更快的收敛。更准确地说，牛顿-拉夫森方法的收敛是二次的，而梯度下降方法的收敛仅是线性的。
www.buymeacoffee.com/guillaumes0
总结来说，本文探讨了牛顿-拉夫森方法这一基本的数值优化技术，重点是其几何解释。方法的优雅和计算能力得到了突出，文章旨在以直观的方式阐明其机制。
牛顿-拉夫森方法被描述为一种求根的迭代过程，其目标是近似一个数学方程的根。核心原则进行了讨论，强调了线性化和切线交点作为关键概念。
文章随后区分了牛顿-拉夫森方法的两个主要应用：求根和优化。求根被解释为寻找函数与 x 轴的交点，而优化被描述为寻找函数导数为零的极值点。方法在优化中的适应性，使用二次多项式来近似函数，被概述了。
应用牛顿-拉夫森优化方法的条件已被提到，包括函数需要是二次可微的要求，以及从一个好的初始近似值开始的重要性。
最后，文章比较了牛顿-拉夫森方法和梯度下降方法在其微分要求、调节和收敛速度方面的区别。文章总结了牛顿-拉夫森方法在优化中的独特优势。
在即将发布的后续文章中，我们将更详细地探讨牛顿-拉夫森方法在实际应用中的应用。
如果你对数学概念感兴趣，并且希望通过 Python 快速学习它们，可以看看我的书：
[## 揭示 70 个数学概念与 Python：通过 Python 探索数学的实用指南…
在…
amzn.to](https://amzn.to/3QaRkXZ?source=post_page-----d9f7acf1b5ae--------------------------------)
优化、牛顿法与利润最大化：第一部分 — 基本优化理论
所有图片由作者提供
学习如何解决和利用牛顿法解决多维优化问题
·
关注 发表在 Towards Data Science ·14 分钟阅读·2023 年 1 月 10 日
–
本文是 3 部分系列中的第一部分。在第一部分中，我们将学习基本的优化理论。然后，在第二部分，我们将扩展这些理论到约束优化问题。最后，在第三部分中，我们将应用所涵盖的优化理论，以及计量经济学和经济理论，来解决一个利润最大化问题。
数学优化是一个极其强大的数学领域，支撑了我们数据科学家在日常工作中隐性或显性使用的许多工具——实际上，几乎所有的机器学习算法都利用优化理论来实现模型收敛。例如，在分类问题中，我们试图通过选择模型的最优参数或权重来最小化对数损失。一般来说，数学优化可以被视为机器学习的主要理论机制。对数学优化的深刻理解是数据科学家工具箱中非常有用的技能——它使数据科学家能够更深入地理解许多当前使用的算法，并且进一步解决各种独特的优化问题。
许多读者可能对梯度下降或相关的优化算法，如随机梯度下降，已经有所了解。然而，本文将更深入地讨论经典的牛顿优化方法，有时称为牛顿-拉夫森方法。我们将从基础数学知识开始，逐步讲解优化理论，到梯度下降，然后深入探讨牛顿方法及其在 Python 中的实现。这将为我们进入第二部分的约束优化和本系列的第三部分中的计量经济学利润最大化问题提供必要的前期准备。
优化基础——一个简单的二次函数
数学优化可以被定义为“确定数学定义问题的最佳解决方案的科学。”[1] 在一些实际例子中，这可以被概念化为：选择参数以最小化机器学习算法的损失函数，选择价格和广告以最大化利润，选择股票以最大化风险调整后的财务回报等等。形式上，任何数学优化问题都可以抽象地表述如下：
(1)
这可以理解为：选择向量x的实值，以最小化目标函数 f(x)（或最大化-f(x))，并满足不等式约束 g(x)和等式约束 h(x)。我们将在本系列的第二部分中讨论如何求解约束优化问题——因为它们使优化问题变得特别复杂。现在，让我们来看一个无约束的单变量示例——考虑以下优化问题：
（2）
在这种情况下，我们想选择使上面的二次函数最小化的 x 值。我们可以采用多种方法——首先，一种简单的方法是对 x 值的大范围进行网格搜索，并选择使f(x)具有最低函数值的 x。然而，随着搜索空间的增加、函数变得更加复杂或维度增加，这种方法可能很快失去计算上的可行性。
如果存在封闭形式的解，我们可以直接使用微积分来求解。也就是说，我们可以通过微积分解析地求解 x 的值。通过取导数（或在高维中，如后面所述的梯度）并将其设置为 0——相对最小值的必要一阶条件——我们可以求解函数的相对极值。然后我们可以取第二导数（或在高维中，如后面所述的 Hessian 矩阵）来确定这个极值是最大值还是最小值。第二导数大于 0（或正定 Hessian 矩阵）——相对最小值的必要二阶条件——意味着是最小值，反之亦然。观察：
（3）
我们可以通过图形验证上述（2）：
请注意，当一个函数存在多个极值点（即多个最小值或最大值）时，必须小心确定哪个是全局极值——我们将在本文中简要讨论这个问题。
上述分析方法可以扩展到更高维度，利用梯度和 Hessian 矩阵——然而，我们不会在高维中求解封闭形式的解——但直觉依然相同。我们仍然会利用迭代方案来解决更高维的问题。我说的迭代方案是什么意思？通常，可能不存在封闭形式（或解析）的解，并且为了存在最大值或最小值，封闭形式的解确实不一定存在。因此，我们需要一种数值方法来解决优化问题。这引导我们到更广泛的迭代方案，包括梯度下降法和牛顿方法。
迭代优化方案
一般来说，迭代优化方案主要分为三类，即 零阶、一阶 和 二阶，它们分别利用函数的零阶、一阶或二阶导数的局部信息。[1] 要使用每种迭代方案，函数 f(x) 必须是相应阶数上连续可微的函数。
零阶迭代方案
零阶迭代方案 与上述的网格搜索紧密相关——简单来说，你在一定范围内搜索可能的 x 值以获得最小的函数值。正如你可能猜到的，这些方法往往比使用高阶方法的计算开销大得多。不用说，它们可以是可靠的并且容易编程。市场上有一些方法可以改进简单的网格搜索，参见[1]了解更多信息；然而，我们将更多关注高阶方案。
一阶迭代方案
一阶迭代方案 是利用目标函数一阶导数的局部信息的迭代方案。最显著的例子是梯度下降方法。对于如上所述的单变量函数，梯度就是一阶导数。将此推广到 n 维度，对于一个函数 f(x)，梯度是一阶偏导数的向量：
(4) 连续可微函数 f(x)
梯度下降从选择一个随机的起点开始，并在 f(x) 的负梯度方向上迭代进行——函数的最陡方向。每次迭代步骤可以表示如下：
(5) 梯度下降迭代方案
其中 γ 是相应的学习率，控制梯度下降算法在每次迭代中“学习”的快慢。学习率过大会导致我们的迭代不受控制地发散。学习率过小则迭代可能需要很长时间才能收敛。此方案会迭代进行，直到达到一个或多个收敛准则，如：
(6) 迭代优化方案的收敛准则
对于某个小的 epsilon 阈值。回到我们的二次例子，将初始猜测设置为 x = 3 和学习率设置为 0.1，步骤如下：
(7)
视觉化如下：
梯度下降和一阶迭代方案在性能上显著可靠。实际上，梯度下降算法主要用于神经网络和机器学习模型中的损失函数优化，许多发展已提高了这些算法的效能。然而，它们仍然仅使用关于函数的有限局部信息（仅一阶导数）。因此，在高维情况下，根据目标函数的性质和学习率，这些方案 1) 可能具有较慢的收敛速度，因为它们保持线性收敛率，并且 2) 可能完全无法收敛。由于这个原因，数据科学家扩大优化工具库是有益的！
二阶迭代方案
如你现在可能已经明白，二阶迭代方案 是利用目标函数的一阶导数和二阶导数的局部信息的迭代方案。最显著的是，我们有牛顿法 (NM)，它使用目标函数的海森矩阵。对于单变量函数，海森矩阵仅仅是二阶导数。类似于梯度，将其推广到 n 维度，海森矩阵是一个 n x n 对称 矩阵，包含一个两次连续可微函数 f(x）的二阶偏导数。
(8) 二次连续可微函数 f(x）的海森矩阵
现在转到导出 NM，首先回忆最小值的一阶必要条件：
(9) x* 处相对最小值的一级必要条件
我们可以使用泰勒级数展开来近似 x*：
(10)
每次迭代的增量 Δ 是对 x* 的一个更好的预期近似。因此，使用 NM 的每次迭代步骤可以表示如下：
(11) 牛顿法迭代方案
回到我们的二次函数示例，将初始猜测设置为 x = 3，步骤如下：
(12)
我们优雅地在第一次迭代时就收敛到最优解。注意，无论方案如何，收敛标准都是相同的。
请注意，所有优化方案都可能陷入局部极值，而不是全局极值（即，考虑具有多个极值（最小值和/或最大值）的高阶多项式——我们可能会陷入一个局部极值，而实际上，另一个极值可能在全球范围内对我们的实际问题更为优化）。已有的方法，并且仍在不断开发，用于处理全局优化，我们将不会深入探讨。你可以利用函数形式的先验知识来设置对结果的预期（即，如果一个严格凸函数有一个临界点，则它必须是全局最小值）。然而，作为一般经验法则，通常明智的做法是对不同的初始值 x 迭代优化方案，然后研究结果的稳定性，通常选择具有最优函数值的结果。
多维示例——罗森布罗克的抛物线谷
现在让我们考虑以下两个变量的优化问题：
(13) 罗森布罗克的抛物线谷
我们将首先通过手动求解上述优化问题，然后在 Python 中进行求解，均使用牛顿法。
通过手动求解
要手动求解，我们需要求解梯度，求解 Hessian，选择我们的初始猜测 Γ = [x,y]，然后迭代将这些信息输入到 NM 算法中，直到收敛为止。首先，求解梯度，我们得到：
(14)
求解 Hessian，我们得到：
(15)
将我们的初始猜测设置为 Γ = [-1.2,1]，我们得到：
(16)
因此，我们成功地求解了目标函数的最优最小值为 Γ* = [1,1]。
通过 Python 求解
我们现在将转向用 Python 求解这个问题，并将其推广到任何函数，使用 SymPy —— 一个用于符号数学的 Python 库。首先，让我们定义罗森布罗克的抛物线谷，并计算该函数的梯度和 Hessian：
import sympy as sm
import numpy as np
# Define symbols & objective function
x, y = sm.symbols('x y')
Gamma = [x,y]
objective = 100*(y-x**2)**2 + (1-x)**2
def get_gradient(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
) -> np.ndarray:
"""
Calculate the gradient of a function.
Args:
function (sm.core.expr.Expr): The function to calculate the gradient of.
symbols (list[sm.core.symbol.Symbol]): The symbols representing the variables in the function.
Returns:
numpy.ndarray: The gradient of the function.
"""
d1 = {}
gradient = np.array([])
for i in symbols:
d1[i] = sm.diff(function, i, 1)
gradient = np.append(gradient, d1[i])
return gradient
def get_hessian(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
) -> np.ndarray:
"""
Calculate the Hessian matrix of a function.
Args:
function (sm.core.expr.Expr): The function for which the Hessian matrix is calculated.
symbols (list[sm.core.symbol.Symbol]): The list of symbols used in the function.
Returns:
numpy.ndarray: The Hessian matrix of the function.
"""
d2 = {}
hessian = np.array([])
for i in symbols:
for j in symbols:
d2[f"{i}{j}"] = sm.diff(function, i, j)
hessian = np.append(hessian, d2[f"{i}{j}"])
hessian = np.array(np.array_split(hessian, len(symbols)))
return hessian
SymPy 允许我们调查方程的符号表示。例如，如果我们调用 objective ，我们将看到相应的输出：
SymPy 的函数符号表示
此外，SymPy 允许我们利用 sm.diff() 命令对相关函数进行求导。如果我们运行定义的函数以获得梯度 get_gradient(objective,Gamma) ，我们得到一个表示梯度的 numpy 数组：
SymPy 求解的梯度
访问特定元素时，我们可以看到符号表示 get_gradient(objective, Gamma)[0]：
SymPy 解出的 df(Γ)/dx
类似地，对于 Hessian 矩阵，我们可以调用 get_hessian(objective, Gamma)：
SymPy 解出的 Hessian 矩阵
访问特定元素 get_hessian(objective,Gamma)[0][1]
SymPy 解出的 df(Γ)/dxdy
可以轻松验证梯度和 Hessian 矩阵与我们手动计算得到的结果是相同的。SymPy 允许对给定符号的指定值评估任何函数。例如，我们可以通过如下调整函数来评估初始猜测处的梯度：
import sympy as sm
import numpy as np
def get_gradient(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float], # Add x0 as argument
) -> np.ndarray:
"""
Calculate the gradient of a function at a given point.
Args:
function (sm.core.expr.Expr): The function to calculate the gradient of.
symbols (list[sm.core.symbol.Symbol]): The symbols representing the variables in the function.
x0 (dict[sm.core.symbol.Symbol, float]): The point at which to calculate the gradient.
Returns:
numpy.ndarray: The gradient of the function at the given point.
"""
d1 = {}
gradient = np.array([])
for i in symbols:
d1[i] = sm.diff(function, i, 1).evalf(subs=x0) # add evalf method
gradient = np.append(gradient, d1[i])
return gradient.astype(np.float64) # Change data type to float
def get_hessian(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
) -> np.ndarray:
"""
Calculate the Hessian matrix of a function at a given point.
Args:
function (sm.core.expr.Expr): The function for which the Hessian matrix is calculated.
symbols (list[sm.core.symbol.Symbol]): The list of symbols used in the function.
x0 (dict[sm.core.symbol.Symbol, float]): The point at which the Hessian matrix is evaluated.
Returns:
numpy.ndarray: The Hessian matrix of the function at the given point.
"""
d2 = {}
hessian = np.array([])
for i in symbols:
for j in symbols:
d2[f"{i}{j}"] = sm.diff(function, i, j).evalf(subs=x0)
hessian = np.append(hessian, d2[f"{i}{j}"])
hessian = np.array(np.array_split(hessian, len(symbols)))
return hessian.astype(np.float64)
我们现在可以通过调用 get_gradient(objective, Gamma, {x:-1.2,y:1.0}) 来计算给定起始点的梯度：
初始点处的梯度
类似地，对于 Hessian 矩阵 get_hessian(objective, Gamma, {x:-1.2,y:1.0})：
初始点处的 Hessian 矩阵
同样，我们可以通过以上手动计算的工作来验证这些值是否正确。现在我们拥有了编写牛顿法所需的所有要素（梯度下降的代码在本文末尾也提供）：
import sympy as sm
import numpy as np
def newton_method(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
iterations: int = 100,
) -> dict[sm.core.symbol.Symbol, float] or None:
"""
Perform Newton's method to find the solution to the optimization problem.
Args:
function (sm.core.expr.Expr): The objective function to be optimized.
symbols (list[sm.core.symbol.Symbol]): The symbols used in the objective function.
x0 (dict[sm.core.symbol.Symbol, float]): The initial values for the symbols.
iterations (int, optional): The maximum number of iterations. Defaults to 100.
Returns:
dict[sm.core.symbol.Symbol, float] or None: The solution to the optimization problem, or None if no solution is found.
"""
x_star = {}
x_star[0] = np.array(list(x0.values()))
print(f"Starting Values: {x_star[0]}")
for i in range(iterations):
gradient = get_gradient(function, symbols, dict(zip(x0.keys(), x_star[i])))
hessian = get_hessian(function, symbols, dict(zip(x0.keys(), x_star[i])))
x_star[i + 1] = x_star[i].T - np.linalg.inv(hessian) @ gradient.T
if np.linalg.norm(x_star[i + 1] - x_star[i]) < 10e-5:
solution = dict(zip(x0.keys(), x_star[i + 1]))
print(f"\nConvergence Achieved ({i+1} iterations): Solution = {solution}")
break
else:
solution = None
print(f"Step {i+1}: {x_star[i+1]}")
return solution
我们现在可以通过 newton_method(objective,Gamma,{x:-1.2,y:1}) 来运行代码：
结论
就这样！如果你已经阅读到这一步，你现在对如何思考和抽象地制定无约束数学优化问题有了扎实的理解，包括基本的分析方法和更复杂的迭代方法。显然，我们在迭代方案中可以融入的信息越多（即更高阶的导数），收敛速度就越高效。请注意，我们只是触及了数学优化复杂世界的表面。 尽管如此，我们今天讨论的工具在实践中绝对可以使用，并可以扩展到更高维的优化问题。
请关注 第二部分 的系列文章，我们将在其中扩展我们在这里学到的内容，解决有约束的优化问题——这是无约束优化的一个极其实用的扩展。实际上，大多数实际优化问题都会对选择变量有某种形式的约束。然后我们将转向本系列的第三部分，在其中我们将应用学到的优化理论和额外的计量经济学与经济理论来解决一个简单的利润最大化问题。希望你和我一样喜欢阅读这篇文章！
奖励——牛顿法的陷阱
尽管牛顿法具有吸引力，但它也有自身的陷阱。尤其是，存在两个主要陷阱——1）即使在选择接近解的起始点时，NM 并不总是收敛；2）NM 在每一步都需要计算 Hessian 矩阵，这在高维情况下计算开销非常大。对于陷阱 #1），一种相应的解决方案是改进牛顿法（MNM），可以粗略地认为它是梯度下降法，其中搜索方向由牛顿步长Δ给出。对于陷阱 #2），已经提出了准牛顿方法，如 DFP 或 BFGS，这些方法通过近似每一步使用的逆 Hessian 矩阵来减轻计算负担。有关更多信息，请参见[1]。
梯度下降函数
import sympy as sm
import numpy as np
def gradient_descent(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
learning_rate: float = 0.1,
iterations: int = 100,
) -> dict[sm.core.symbol.Symbol, float] or None:
"""
Performs gradient descent optimization to find the minimum of a given function.
Args:
function (sm.core.expr.Expr): The function to be optimized.
symbols (list[sm.core.symbol.Symbol]): The symbols used in the function.
x0 (dict[sm.core.symbol.Symbol, float]): The initial values for the symbols.
learning_rate (float, optional): The learning rate for the optimization. Defaults to 0.1.
iterations (int, optional): The maximum number of iterations. Defaults to 100.
Returns:
dict[sm.core.symbol.Symbol, float] or None: The solution found by the optimization, or None if no solution is found.
"""
x_star = {}
x_star[0] = np.array(list(x0.values()))
print(f"Starting Values: {x_star[0]}")
for i in range(iterations):
gradient = get_gradient(function, symbols, dict(zip(x0.keys(), x_star[i])))
x_star[i + 1] = x_star[i].T - learning_rate * gradient.T
if np.linalg.norm(x_star[i + 1] - x_star[i]) < 10e-5:
solution = dict(zip(x0.keys(), x_star[i + 1]))
print(f"\nConvergence Achieved ({i+1} iterations): Solution = {solution}")
break
else:
solution = None
print(f"Step {i+1}: {x_star[i+1]}")
return solution
资源
[1] Snyman, J. A., & Wilke, D. N. (2019). 实用数学优化：基本优化理论与基于梯度的算法（第 2 版）。Springer。
[2] en.wikipedia.org/wiki/Gradient_descent
通过此 GitHub Repo 访问所有代码: github.com/jakepenzak/Blog-Posts
感谢你阅读我的帖子！我在 Medium 上的帖子旨在探讨利用 计量经济学 和 统计/机器学习 技术的现实世界和理论应用。此外，我还希望通过理论和模拟提供某些方法论的理论基础。最重要的是，我写作是为了学习！我希望让复杂的话题对所有人稍微更易于理解。如果你喜欢这篇文章，请考虑 在 Medium 上关注我！
优化、牛顿法与利润最大化：第二部分——约束优化理论
作者提供的所有图片
了解如何扩展牛顿法并解决约束优化问题
·
关注 发表在 Towards Data Science ·13 分钟阅读·2023 年 2 月 2 日
–
这篇文章是一个三部分系列中的第二部分。在第一部分中，我们研究了基本的优化理论。现在，在第二部分中，我们将把这一理论扩展到受限优化问题。最后，在第三部分中，我们将应用所涵盖的优化理论，以及计量经济学和经济理论，来解决一个利润最大化问题。
考虑以下问题：你想确定在特定金融工具上投资多少资金以最大化你的投资回报。然而，单纯最大化投资回报的问题过于宽泛和简单。由于其简单性，解决方案就是将所有资金投入到回报潜力最高的金融工具中。显然，这不是一个好的投资策略；那么，我们如何改进呢？通过对投资决策施加约束，我们的选择变量。 例如，我们可以指定约束条件，比如 1) 限制我们愿意承担的金融风险量（见现代投资组合理论），或者 2) 指定我们投资组合中每类金融工具（如股票、债券、衍生品等）的分配比例——可能性无穷无尽。注意，当我们添加约束时，这个问题变得显著更具可处理性。尽管这是一个简单的例子，但它有助于捕捉受限优化的一个基本动机：
受限优化的本质在于为无约束优化问题提供一种适用性和解决复杂现实世界问题的能力。
受限优化被定义为“在变量受到约束的情况下，对目标函数进行优化的过程。”[1] 添加对变量的约束将一个无约束的、或许是难以处理的优化问题转化为一个有助于建模和解决现实世界问题的问题。然而，添加约束可能将一个简单的优化问题转变为一个不再是微不足道的问题。在这篇文章中，我们将深入探讨一些可以添加到我们工具箱中的技术，以扩展在第一部分中学习的无约束优化理论，从而解决受限优化问题。
在 第一部分中，我们介绍了基本的优化理论——包括 1) 解析设置和解决一个简单的单变量优化问题，2) 迭代优化方案——即梯度下降法和牛顿法，以及 3) 手动和使用 Python 实现牛顿法用于多维优化问题。本文旨在使那些已经熟悉 第一部分 中内容的读者能够轻松理解。
约束优化基础（& 第一部分 回顾）
一个数学优化问题可以抽象地表示如下：
(1)
在这里，我们选择 x 的实际值，以最小化 目标函数 f(x)（或最大化 -f(x）），同时满足 不等式约束 g(x) 和 等式约束 h(x)。在 第一部分 中，我们讨论了如何在没有 g(x) 和 h(x) 的情况下解决这些问题，现在我们将这些约束重新引入到我们的优化问题中。首先，让我们简明扼要地回顾如何实现牛顿法用于无约束问题。
回顾一下，我们可以使用泰勒级数展开来近似最小值的一阶必要条件：
(2)
其中 H(x) 和 ∇f(x) 分别表示 f(x) 的 Hessian 矩阵和梯度。每次迭代增加的 delta, Δ, 是对最优值 x* 的预期更好近似。因此，每次使用牛顿法的迭代步骤可以表示如下：
(3) 牛顿法迭代方案
我们执行这个方案直到在以下一个或多个标准上达到收敛：
(4) 迭代优化方案的收敛标准
将其转化为 Python 代码，我们使用 SymPy —— 一个用于符号数学的 Python 库 —— 并创建可泛化的函数来计算梯度、计算 Hessian 矩阵，并实现牛顿法用于 n 维函数：
import sympy as sm
import numpy as np
def get_gradient(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
) -> np.ndarray:
"""
Calculate the gradient of a function at a given point.
Args:
function (sm.core.expr.Expr): The function to calculate the gradient of.
symbols (list[sm.core.symbol.Symbol]): The symbols representing the variables in the function.
x0 (dict[sm.core.symbol.Symbol, float]): The point at which to calculate the gradient.
Returns:
numpy.ndarray: The gradient of the function at the given point.
"""
d1 = {}
gradient = np.array([])
for i in symbols:
d1[i] = sm.diff(function, i, 1).evalf(subs=x0)
gradient = np.append(gradient, d1[i])
return gradient.astype(np.float64)
def get_hessian(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
) -> np.ndarray:
"""
Calculate the Hessian matrix of a function at a given point.
Args:
function (sm.core.expr.Expr): The function for which the Hessian matrix is calculated.
symbols (list[sm.core.symbol.Symbol]): The list of symbols used in the function.
x0 (dict[sm.core.symbol.Symbol, float]): The point at which the Hessian matrix is evaluated.
Returns:
numpy.ndarray: The Hessian matrix of the function at the given point.
"""
d2 = {}
hessian = np.array([])
for i in symbols:
for j in symbols:
d2[f"{i}{j}"] = sm.diff(function, i, j).evalf(subs=x0)
hessian = np.append(hessian, d2[f"{i}{j}"])
hessian = np.array(np.array_split(hessian, len(symbols)))
return hessian.astype(np.float64)
def newton_method(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
iterations: int = 100,
) -> dict[sm.core.symbol.Symbol, float] or None:
"""
Perform Newton's method to find the solution to the optimization problem.
Args:
function (sm.core.expr.Expr): The objective function to be optimized.
symbols (list[sm.core.symbol.Symbol]): The symbols used in the objective function.
x0 (dict[sm.core.symbol.Symbol, float]): The initial values for the symbols.
iterations (int, optional): The maximum number of iterations. Defaults to 100.
Returns:
dict[sm.core.symbol.Symbol, float] or None: The solution to the optimization problem, or None if no solution is found.
"""
x_star = {}
x_star[0] = np.array(list(x0.values()))
# x = [] ## Return x for visual!
print(f"Starting Values: {x_star[0]}")
for i in range(iterations):
# x.append(dict(zip(x0.keys(),x_star[i]))) ## Return x for visual!
gradient = get_gradient(function, symbols, dict(zip(x0.keys(), x_star[i])))
hessian = get_hessian(function, symbols, dict(zip(x0.keys(), x_star[i])))
x_star[i + 1] = x_star[i].T - np.linalg.inv(hessian) @ gradient.T
if np.linalg.norm(x_star[i + 1] - x_star[i]) < 10e-5:
solution = dict(zip(x0.keys(), x_star[i + 1]))
print(f"\nConvergence Achieved ({i+1} iterations): Solution = {solution}")
break
else:
solution = None
print(f"Step {i+1}: {x_star[i+1]}")
return solution
为了解决无约束优化问题，我们可以运行以下代码：
import sympy as sm
# Define Symbols
x, y = sm.symbols('x y')
Gamma = [x,y]
# Define Objective Function (Rosenbrock's Parabolic Valley)
objective = 100*(y-x**2)**2 + (1-x)**2
# Specify starting values
Gamma0 = {x:-1.2,y:1}
# Call function
newton_method(objective, Gamma, Gamma0)
及其对应的输出：
如果上述所有材料感觉非常陌生，那么我建议查看part 1，它将更深入地探讨上述内容并帮助你跟上进度！话不多说，让我们深入实施优化问题中的约束。
注意：所有以下约束优化技术都可以且应该在适用时与梯度下降算法结合使用！
求解带约束的优化问题
如我们上面讨论的，目标函数可能有两种约束——等式约束和不等式约束。注意，对于每种类型的约束，都有不同的方法，具有不同的优缺点。有关不同方法的进一步讨论，请参见[2]。不过，我们将重点关注两种方法，一种用于等式约束，另一种用于不等式约束，这些方法在性能上可靠，对新手易于理解，并且可以很容易地集成到一个有机的问题中。
等式约束 — 拉格朗日
首先，我们将处理带有等式约束的优化问题。即，具有以下形式的优化问题：
(5)
假设我们正在处理 Rosenbrock 的抛物谷，如part 1中所述，但现在添加了等式约束 x² - y = 2：
(6) 带有等式约束的 Rosenbrock 的抛物谷 （问题 1）
注意，为了简化和一致性，等式约束应写成等于零的形式。现在我们的优化问题看起来像：
Rosenbrock 的抛物谷（紫黄色色图）和等式约束曲线（黑色）
在这里，可行区域的最优值位于等式约束曲线与我们上面的目标函数的交点之一。
约瑟夫-路易·拉格朗日开发了一种方法，将等式约束直接纳入目标函数中——创建拉格朗日函数——以便可以继续应用传统方法使用一阶和二阶导数。[2][3] 形式上，拉格朗日函数具有以下形式：
(7) 拉格朗日函数的正式定义
其中 f(x) 和 h(x) 分别是目标函数和等式约束。Λ 是与每个等式约束 j 对应的 拉格朗日乘子。拉格朗日乘子被视为拉格朗日函数中的新选择变量。正好，x* 作为等式约束问题的最小值的必要条件是 x* 对应于拉格朗日函数的驻点 (x*, Λ*). 即，
(8) 拉格朗日一阶条件
对于上述示例——等式约束的 Rosenbrock 抛物线谷（公式 1）——我们可以将拉格朗日函数写为：
(9)
然后我们可以使用牛顿法求解这个拉格朗日函数，但现在需要将拉格朗日乘子作为附加选择变量。
import sympy as sm
x, y, λ = sm.symbols('x y λ')
Langrangian_objective = 100*(y-x**2)**2 + (1-x)**2 + λ*(x**2-y-2)
Gamma = [x,y,λ]
Gamma0 = {x:-1.2,y:1,λ:1}
newton_method(Langrangian_objective,Gamma,Gamma0)
对应的输出为：
可以很容易验证解满足我们的等式约束。就这样！这还不算太难，对吧？这种方法可以扩展以添加任何数量的等式约束——只需添加另一个拉格朗日乘子。现在我们继续讨论如何纳入不等式约束。
不等式约束—对数障碍函数
现在我们将处理带有不等式约束的优化问题。即，具有如下形式的优化问题：
(10)
再次假设，我们在处理 Rosenbrock 的抛物线谷，但现在有不等式约束 x ≤ 0 和 y ≥ 3：
(11) 带有不等式约束的 Rosenbrock 抛物线谷 （问题 2）
现在我们的优化问题变成了：
Rosenbrock 的抛物线谷（紫黄色色彩图）和不等式约束平面（黑色）
最优值的可行区域位于由红星标记的约束所界定的象限中。
由于这些约束没有严格的等式，我们无法将它们直接纳入目标函数。然而，我们可以动脑筋——我们可以做的是增强目标函数，在目标函数中加入一个“障碍”，对接近不等式约束边界的解值进行惩罚。这些方法被称为“内点法”或“障碍法”[4][5]。像拉格朗日函数一样，我们可以通过引入障碍函数（在我们这个案例中是对数障碍函数）将原来的约束优化问题转化为无约束优化问题，从而创建障碍函数。形式上，对数障碍函数的特点是：
(12) 对数障碍函数的正式定义
其中ρ是一个小的正标量——称为障碍参数。随着ρ → 0，障碍函数B(x,ρ)的解应收敛到我们原始约束优化函数的解。注意，c(x)表示，根据我们如何制定不等式约束（大于或小于零），将决定我们使用该约束的负值或正值。我们知道 y=log(x)在 x ≤ 0 时未定义，因此我们需要将约束制定为始终≥0。
你可能会问，对数障碍法究竟是如何工作的？首先，在使用障碍法时，我们必须选择位于可行区域的起始值。随着最优值接近由约束定义的“障碍”，该方法依赖于对数函数在值接近零时趋向负无穷的特性，从而惩罚目标函数值。随着ρ → 0，惩罚减小（见下图），我们逐渐收敛到解。然而，需要从足够大的ρ开始，以确保惩罚足够大，防止“跳出”障碍。因此，该算法比牛顿法多了一个额外的循环——即，我们选择一个起始值ρ，使用牛顿法优化障碍函数，然后通过缓慢减小ρ（ρ → 0）来更新ρ，直到收敛。
不同值的ρ的对数障碍
回到我们之前的例子——不等式约束的罗森布罗克抛物谷（公式 2）——我们可以将障碍函数写为：
(13)
记住 log(a) + log(b) = log(ab)，以及我们的一个约束 x ≤ 0 → -x ≥ 0。我们必须更新我们的代码以适应障碍法算法：
import sympy as sm
import numpy as np
def constrained_newton_method(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
iterations: int = 100,
) -> dict[sm.core.symbol.Symbol, float] or None:
"""
Performs constrained Newton's method to find the optimal solution of a function subject to constraints.
Parameters:
function (sm.core.expr.Expr): The function to optimize.
symbols (list[sm.core.symbol.Symbol]): The symbols used in the function.
x0 (dict[sm.core.symbol.Symbol, float]): The initial values for the symbols.
iterations (int, optional): The maximum number of iterations. Defaults to 100.
Returns:
dict[sm.core.symbol.Symbol, float] or None: The optimal solution if convergence is achieved, otherwise None.
"""
x_star = {}
x_star[0] = np.array(list(x0.values())[:-1])
optimal_solutions = []
optimal_solutions.append(dict(zip(list(x0.keys())[:-1], x_star[0])))
for step in range(iterations):
# Evaluate function at rho value
if step == 0: # starting rho
rho_sub = list(x0.values())[-1]
rho_sub_values = {list(x0.keys())[-1]: rho_sub}
function_eval = function.evalf(subs=rho_sub_values)
print(f"Step {step} w/ {rho_sub_values}") # Barrier method step
print(f"Starting Values: {x_star[0]}")
# Newton's Method
for i in range(iterations):
gradient = get_gradient(
function_eval, symbols[:-1], dict(zip(list(x0.keys())[:-1], x_star[i]))
)
hessian = get_hessian(
function_eval, symbols[:-1], dict(zip(list(x0.keys())[:-1], x_star[i]))
)
x_star[i + 1] = x_star[i].T - np.linalg.inv(hessian) @ gradient.T
if np.linalg.norm(x_star[i + 1] - x_star[i]) < 10e-5:
solution = dict(zip(list(x0.keys())[:-1], x_star[i + 1]))
print(
f"Convergence Achieved ({i+1} iterations): Solution = {solution}\n"
)
break
# Record optimal solution & previous optimal solution for each barrier method iteration
optimal_solution = x_star[i + 1]
previous_optimal_solution = list(optimal_solutions[step - 1].values())
optimal_solutions.append(dict(zip(list(x0.keys())[:-1], optimal_solution)))
# Check for overall convergence
if np.linalg.norm(optimal_solution - previous_optimal_solution) < 10e-5:
print(
f"\n Overall Convergence Achieved ({step} steps): Solution = {optimal_solutions[step]}\n"
)
overall_solution = optimal_solutions[step]
break
else:
overall_solution = None
# Set new starting point
x_star = {}
x_star[0] = optimal_solution
# Update rho
rho_sub = 0.9 * rho_sub
return overall_solution
我们现在可以用上述代码求解障碍函数（注意：确保起始值在不等式约束的可行范围内，如果跳出不等式约束，可能需要增加ρ的起始值）：
import sympy as sm
x, y, ρ = sm.symbols('x y ρ')
Barrier_objective = 100*(y-x**2)**2 + (1-x)**2 - ρ*sm.log((-x)*(y-3))
Gamma = [x,y,ρ] # Function requires last symbol to be ρ!
Gamma0 = {x:-15,y:15,ρ:10}
constrained_newton_method(Barrier_objective,Gamma,Gamma0)
其对应的输出为：
很明显，解满足指定的不等式约束。就是这样。我们现在已经处理了优化问题中的不等式约束。最后，让我们将一切整合起来，继续处理具有混合约束的约束优化问题——这只是我们上述工作组合的结果。
汇总
现在，让我们通过结合上述的等式和不等式约束来解决我们的优化问题。也就是说，我们想要解决一个形式的优化问题：
(14)
我们只需要将拉格朗日函数和障碍函数结合成一个函数。因此，我们可以创建一个通用函数，称之为 O，用于处理具有等式和不等式约束的优化问题：
(15) 可推广的约束优化问题的函数
其中，如前所述，Λ 是拉格朗日乘子向量，ρ 是障碍参数。因此，结合我们上面的受限（公式 6）和无约束问题（公式 11），我们可以将我们的混合受限优化问题表述如下：
(16)
在 Python 中，
import sympy as sm
x, y, λ, ρ = sm.symbols('x y λ ρ')
combined_objective = 100*(y-x**2)**2 + (1-x)**2 + λ*(x**2-y-2) - ρ*sm.log((-x)*(y-3))
Gamma = [x,y,λ,ρ] # Function requires last symbol to be ρ!
Gamma0 = {x:-15,y:15,λ:0,ρ:10}
constrained_newton_method(combined_objective,Gamma,Gamma0)
以及相应的输出：
我们可以验证这个解确实符合约束条件。具体来说，x ≤ 0，y ≥ 3，& x² - y = 2。令人满意，不是吗？
结论
哎呀。深呼吸一下——你值得拥有。希望到目前为止，你对将约束条件融入优化问题的技术有了更好的理解。我们仍然只是触及了数学优化中各种工具和技术的表面。
敬请关注系列的第三部分，即最后一部分，我们将应用到目前为止学到的优化材料，并结合计量经济学和经济理论来解决利润最大化问题。我的目标是第三部分将总结我们所涵盖的内容，并展示一个实际应用案例。像往常一样，我希望你能像我写作时一样享受阅读这篇文章！
资源
[1] en.wikipedia.org/wiki/Constrained_optimization
[2] Snyman, J. A., & Wilke, D. N. (2019). 实用数学优化：基本优化理论与基于梯度的算法（第 2 版）。Springer。
[3] en.wikipedia.org/wiki/Lagrange_multiplier
[4] en.wikipedia.org/wiki/Interior-point_method
[5] en.wikipedia.org/wiki/Barrier_function
通过这个 GitHub 仓库访问所有代码： github.com/jakepenzak/Blog-Posts
感谢你阅读我的文章！我在 Medium 上的文章旨在探索利用 计量经济学 和 统计/机器学习 技术的现实世界和理论应用。此外，我还希望通过理论和模拟提供各种方法论的理论基础。最重要的是，我写作是为了学习！我希望能让复杂的话题对所有人更易于理解。如果你喜欢这篇文章，请考虑 在 Medium 上关注我！
优化、牛顿法与利润最大化：第三部分 — 应用利润最大化
所有图片由作者提供
学习如何应用优化和计量经济学技术来解决实际利润最大化问题
·
关注 发表在 Towards Data Science ·19 分钟阅读·2023 年 3 月 1 日
–
本文是第三篇，也是最后一篇，在这个三部分系列中。在第一部分中，我们研究了基本优化理论。然后，在第二部分中，我们将这一理论扩展到约束优化问题。现在，在第三部分中，我们将应用所涉及的优化理论，以及计量经济学和经济学理论，来解决利润最大化问题。
假设，作为公司的一名数据科学家，你的任务是估算最佳的资金分配，以最大化某产品线的整体利润。此外，假设你在这些分配决策上有一些约束条件，例如必须分配的总花费上限和/或必须在某些渠道上花费的最低金额。在本文中，我们将结合第一部分和第二部分中涉及的优化理论，以及额外的经济学和计量经济学理论，来处理这种理论上的利润最大化问题——我们将在本文中进一步详细阐述。
本文的目标是将我们迄今为止学到的知识结合起来，我希望能激励和启发读者如何将这些技术应用到实际环境中。本文并非针对所讨论问题的全面解决方案，因为细微差别和特性当然会使理论示例复杂化。此外，许多涉及的技术在 Python 中通过像pyomo、SciPy等包有更优化的实现。尽管如此，我希望能提供一个构建应用优化问题的强大框架。让我们深入探讨吧！
在第一部分中，我们涵盖了基本的优化理论 — 包括 1) 设置并通过分析解决一个简单的单变量优化问题，2) 迭代优化方案 — 即，梯度下降和牛顿法，和 3) 通过手工和 python 实现牛顿法解决多维优化问题。在第二部分中，我们涵盖了有约束优化理论 — 包括 1)将等式约束和 2)将不等式约束纳入我们的优化问题中，并通过牛顿法解决它们。本文旨在为那些已经熟悉第一部分和第二部分中涵盖的内容的人士提供帮助。*
优化理论（第一部分 & 第二部分 总结）
一个数学优化问题可以抽象地表述如下：
(1)
我们选择使向量x的实值最小化目标函数f(x)（或最大化-f*(x)），受不等式约束g(x)和等式约束h(x)约束。在第二部分中，我们讨论了如何将这些约束直接纳入我们的优化问题中。值得注意的是，使用拉格朗日乘子和对数障碍函数，我们可以构建一个新的目标函数O(x, Λ,** ρ):
(2) 有约束优化问题的泛化函数
其中Λ是与每个等式约束 h(x)相关的拉格朗日乘子向量，而ρ是与所有不等式约束 g(x)相关的障碍参数。然后我们可以通过选择一个初始值ρ（注意，目标函数的大功能值会要求更大的初始值ρ以调整惩罚），使用牛顿法迭代优化新目标函数，接着通过逐渐减小ρ（ρ → 0）来更新ρ，并重复直到收敛 — 牛顿法迭代方案如下：
(3) 牛顿法迭代方案
其中H(x)和∇f(x)分别表示我们的目标函数 O(x, Λ,** ρ)的海森矩阵和梯度。当我们达到以下一个或多个标准的收敛时，即表示收敛：
(4) 迭代优化方案的收敛标准
在 Python 中，利用SymPy，我们有 4 个函数。一个函数用于获取我们的 SymPy 函数的梯度，一个函数用于获取我们的 SymPy 函数的海森矩阵，一个函数通过牛顿法解决无约束优化问题，另一个函数通过牛顿法解决根据广义方程（2）的有约束优化问题（有关该材料的更详细覆盖，请参见第一部分和第二部分）。有关这些函数，请参见附录。
要解决一个有约束的优化问题，我们可以运行以下代码（确保起始值在不等式约束的可行范围内！）：
import sympy as sm
x, y, λ, ρ = sm.symbols('x y λ ρ')
# f(x): 100*(y-x**2)**2 + (1-x)**2
# h(x): x**2 - y = 2
# g_1(x): x <= 0
# g_2(x) y >= 3
objective = 100*(y-x**2)**2 + (1-x)**2 + λ*(x**2-y-2) - ρ*sm.log((y-3)*(-x))
symbols = [x,y,λ,ρ] # Function requires last symbol to be ρ (barrier parameter)
x0 = {x:-15,y:15,λ:0,ρ:10}
constrained_newton_method(objective,symbols,x0)
对应的输出为：
如果上述材料感觉陌生或需要更详细的回顾，建议查看本系列的第一部分和第二部分，这些将提供更深入的材料概述。对于本文的其余部分，我们将首先讨论基本的利润最大化和计量经济学理论，然后进入理论示例的解决方案。
应用利润最大化
假设我们有一个$100,000 的广告预算，并且必须全部花费。我们的任务是选择最优的预算分配到两种广告渠道（数字广告和电视广告），以最大化某一产品线的整体利润。此外，假设我们必须至少分配$20,000 到电视广告和$10,000 到数字广告。
理论公式
现在我们将数学地制定我们寻求解决的利润最大化问题：
(5) 链接
*π(⋅)表示利润函数，δ表示数字广告支出，τ表示电视广告支出，(⋅)是额外变量的占位符。请注意，我们在最小化π(⋅)的负值，这等同于最大化π(*⋅)。利润函数定义如下：
(6) 链接
p表示价格，q(δ, τ, ⋅)表示需求量函数，C*(q(⋅), δ, τ)表示成本函数，直观上，它是数量的直接函数（如果我们生产更多，成本会增加）以及广告支出的多少。成本函数还可以接受额外输入，但为了演示，我们将其保持为数量和广告费用的函数。请注意，我们选择的δ和τ*通过其对需求量和成本函数的影响直接影响利润函数。为了增加优化问题的可操作性，我们需要使用计量经济学技术来估计我们的需求量函数。一旦我们指定了成本函数并估计了需求量函数，我们可以按如下方式解决优化问题：
(7) 链接
其中 q_hat 是我们对需求量的经济计量模型的估计。在我们阐述需求量模型的经济计量规格之前，有必要讨论一个关于优化问题所需假设的重要说明，以证明其可操作性。我们必须获得数字广告和电视广告对需求量的 因果估计。用经济学术语来说，数字广告和电视广告在经济计量模型中需要是 外生 的。也就是说，它们与模型中的误差不相关。外生性可以通过两种方式实现：1）我们有正确的经济计量模型结构规范，以衡量数字广告和电视广告对需求量的影响（即，我们包括所有与需求量以及数字和电视广告支出相关的变量）；或 2）我们有随机变化的数字和电视广告支出（这可以通过随时间随机变化的支出来观察需求的响应）。
直观上，外生性是必需的，因为需要捕捉广告支出变化的因果影响——也就是说，如果我们改变广告支出的数值，平均而言会发生什么。如果我们估计的效果 不是 因果的，那么我们对广告支出的变化 不会 对应于需求量的真实变化。注意，模型不需要为需求量做出最佳预测，而是需要准确捕捉因果关系。请参阅我的 前一篇文章，其中更深入地讨论了外生性假设。
现在假设我们为按时间 t 索引的需求量指定以下经济计量模型：
(8) 链接
其中 β 和 γ 分别是数字广告支出自然对数的影响估计值 δ 和电视广告支出的影响估计值 τ。此外，α 是我们的截距，ϕ1 和 ϕ2 是需求量的 自回归 组件的估计值，S 表示季节性，X 是所有相关协变量和滞后协变量的集合及其系数估计矩阵 Ω，ϵ 是误差项。此外，假设数字和电视广告在 X、S 和我们模型中的自回归组件 条件下 满足外生性假设。也就是说，
(9) 链接
你可能会问，为什么使用数字和电视广告支出的自然对数？这并不是在此背景下的必需或决定性的选择，但我旨在展示变量变换如何捕捉我们选择变量与感兴趣结果之间的关系。以我们的案例为例，假设我们假设广告支出的影响最初急剧增加，但逐渐趋于平稳（即收益递减法则）。这正是对数变换允许我们建模的内容。请观察：
注意，成本函数的形式通常可以在内部更容易得知。因此，我们也来指定我们的成本函数的形式：
(10) 链接
在这里，我们可以看到每生产一个单位都有一个成本ζ，并且随着生产量的增加，这个成本会得到折扣（可以理解为对大合同的折扣或规模经济）。我们还简单地将数字广告支出和电视广告支出合并到我们的总成本中。
现在我们已经为我们的计量经济学利润最大化问题建立了理论基础，让我们来模拟一些数据并将其带入 Python！
数据模拟（可选）
请注意，本节可以跳过而不会丢失主要内容。
首先，我们来模拟 10 年期间的每月需求量数据，其中包括以下变量：
import pandas as pd
df = pd.DataFrame()
## Digital Advertising - ln(δ)
df['log_digital_advertising'] = np.log(np.random.normal(loc=50000,scale=15000,size=120).round())
## Television Advertising - ln(τ)
df['log_television_advertising'] = np.log(np.random.normal(loc=50000,scale=15000,size=120).round())
## Matrix X of covariates
# Lag Digital Advertising
df['log_digital_advertising_lag1'] = df['log_digital_advertising'].shift(1)
df['log_digital_advertising_lag2'] = df['log_digital_advertising'].shift(2)
# Lag Television Advertising
df['log_television_advertising_lag1'] = df['log_television_advertising'].shift(1)
df['log_television_advertising_lag2'] = df['log_television_advertising'].shift(2)
# Price
df['price'] = np.random.normal(loc=180,scale=15,size=120).round()
df['price_lag1'] = df['price'].shift(1)
df['price_lag2'] = df['price'].shift(2)
# Competitor Price
df['comp_price'] = np.random.normal(loc=120,scale=15,size=120).round()
df['comp_price_lag1'] = df['comp_price'].shift(1)
df['comp_price_lag2'] = df['comp_price'].shift(2)
# Seasonality
from itertools import cycle
months = cycle(['Jan','Feb','Mar','Apr','May','June','July','Aug','Sep','Oct','Nov','Dec'])
df['months'] = [next(months) for m in range(len(df))]
one_hot = pd.get_dummies(df['months'], dtype=int)
one_hot = one_hot[['Jan','Feb','Mar','Apr','May','June','July','Aug','Sep','Oct','Nov','Dec']]
df = df.join(one_hot).drop('months',axis=1)
## Constant
df['constant'] = 1
# Drop NaN (Two lags)
df = df.dropna()
注意，我们包括滞后变量，因为今天的需求量很可能是许多变量滞后值的函数。我们还通过引入每月的虚拟变量来控制季节性效应（这只是将季节性纳入模型的多种方式之一）。然后我们指定与每个变量相关的参数（请注意，这些参数按数据框的列顺序指定！）：
params = np.array(
[10_000, # β
5_000, # γ
2_000, # Ω
1_000, # Ω
3_000, # Ω
1_000, # Ω
-1_000, # Ω
-500, # Ω
-100, # Ω
500, # Ω
300, # Ω
100, # Ω
25_000, # S
15_000, # S
15_000, # S
10_000, # S
10_000, # S
10_000, # S
15_000, # S
15_000, # S
25_000, # S
35_000, # S
35_000, # S
40_000, # S
50_000 # α
])
然后我们可以通过运行quantity_demanded = np.array(df) @ params来模拟我们的计量经济学规范（方程 8）中的需求量。然而，请注意我们缺少自回归成分，因此我们还希望需求量遵循上述的自回归过程。也就是说，需求量也是其自身滞后值的函数。我们在这里包括 2 个滞后期（AR(2)过程），其系数分别为ϕ1和ϕ2。注意，我们可以通过以下系统用初始条件q0和q-1来模拟：
(11) 链接
def quantity_ar2_process(T, ϕ1, ϕ2, q0, q_1, ϵ, df, params):
Φ = np.identity(T) # The T x T identity matrix
for i in range(T):
if i-1 >= 0:
Φ[i, i-1] = -ϕ1
if i-2 >= 0:
Φ[i, i-2] = -ϕ2
B = np.array(df) @ params + ϵ
B[0] = B[0] + ϕ1 * q0 + ϕ2 * q_1
B[1] = B[1] + ϕ2 * q0
return np.linalg.inv(Φ) @ B
## Quantity Demand AR(2) component process
# Parameters
T = 118 # Time periods less two lags
ϕ1 = 0.3 # Lag 1 coefficient (ϕ1)
ϕ2 = 0.05 # Lag 2 coefficient (ϕ2)
q_1 = 250_000 # Initial Condition q_-1
q0 = 300_000 # Initial Condition q_0
ϵ = np.random.normal(0, 5000, size=T) # Random Error (ϵ)
quantity_demanded_ar = quantity_ar2_process(T,ϕ1,ϕ2,q0,q_1,ϵ,df,params)
# Quantity_demanded target variable
df['quantity_demanded'] = quantity_demanded_ar
# Additional covariates of lagged quantity demanded
df['quantity_demanded_lag1'] = df['quantity_demanded'].shift(1)
df['quantity_demanded_lag2'] = df['quantity_demanded'].shift(2)
计量经济学估计与优化
首先使用我们在方程（2）中的框架，将我们在方程（7）中的约束优化问题转化为一个可以利用我们上面提到的函数constrained_newton_method()来解决的问题：
(12) 链接
如前所述，我们需要估计我们的需求量，q_hat。让我们看看在模拟的 10 年中我们的需求量是什么样的：
我们可以清楚地看到，在年末时出现了一些季节性现象，且我们似乎在处理一个平稳过程（这一切都是由构造决定的）。现在假设我们有以下观察变量：
在方程 8 中，我们的计量经济学规范中，quantity_demanded 是我们的结果 q，log_digital_advertising 是我们的 ln(δ)，log_television_advertising 是我们的 ln(τ)，constant 是我们的 α，quantity_demanded_lag1 和 quantity_demanded_lag2 是我们的自回归组件 q_t-1 和 q_t-2，其余是我们的额外协变量 X，包括季节性 S。
现在，有了这些数据，我们试图估计方程 8 中的计量经济学规范。我们可以使用最小二乘法（OLS）来估计这个结构模型。为此，我们将使用statsmodels。
一个很好的练习是使用我们构建的牛顿法代码来解决线性回归，并将结果与 statsmodels 进行比较。提示：线性回归中的目标是最小化残差平方和。请注意，我们编写的代码绝非解决线性回归的高效方法，但它更多地用于在模型拟合（回归）背景下说明优化理论。相关代码将在文章末尾提供！
注意，我们删除了前两个观察值，因为这些是我们的前两个滞后项，同时我们删除了七月作为参考月份：
import statsmodels.api as stats
## Fit Econometric model using OLS
df = df[2:] # Drop first two lagged values
y = df['quantity_demanded']
X = df.drop(['quantity_demanded','July'],axis=1)
mod = stats.OLS(y,X)
results = mod.fit()
print(results.summary())
现在我们有了我们对需求量的估计计量经济学规范！几点观察：
-
数字广告支出和电视广告支出的对数增加与需求量的增加相关
-
价格的增加与需求量的减少相关（这是预期的行为）
-
我们看到在 9 月到 12 月期间有明显的季节性需求增长，这与我们上面的时间序列一致
-
我们看到需求量的第一个滞后项对现在有预测作用，这支持自回归过程
- 上面的结果可以通过数据模拟部分的数据显示和比较来验证
现在让我们为优化问题指定我们的符号变量（δ、τ、λ 和 ρ），设定时间 t 上的当前变量值，并从数据中获取滞后值。这样我们就有了构建优化问题所需的一切：
# Build Symbolic Functions with all variables in function
δ, τ, λ, ρ = sm.symbols('δ τ λ ρ')
## Values of current variables
price = 180
comp_price = 120
Jan = 1
## Obtain Lagged Values
log_digital_advertising_lag1 = df['log_digital_advertising_lag1'].iloc[-1]
log_digital_advertising_lag2 = df['log_digital_advertising_lag2'].iloc[-2]
log_television_advertising_lag1 = df['log_television_advertising_lag1'].iloc[-1]
log_television_advertising_lag2 = df['log_television_advertising_lag2'].iloc[-2]
price_lag1 = df['price_lag1'].iloc[-1]
price_lag2 = df['price_lag2'].iloc[-2]
comp_price_lag1 = df['comp_price_lag1'].iloc[-1]
comp_price_lag2 = df['comp_price_lag2'].iloc[-2]
quantity_demanded_lag1 = df['quantity_demanded_lag1'].iloc[-1]
quantity_demanded_lag2 = df['quantity_demanded_lag2'].iloc[-2]
variables = [sm.log(δ),
sm.log(τ),
log_digital_advertising_lag1,
log_digital_advertising_lag2,
log_television_advertising_lag1,
log_television_advertising_lag2,
price,
price_lag1,
price_lag2,
comp_price,
comp_price_lag1,
comp_price_lag2,
Jan,0,0,0,0,0,0,0,0,0,0, # All Months less July
1, # Constant
quantity_demanded_lag1,
quantity_demanded_lag2
]
# Quantity Demanded
quantity_demanded = np.array([variables]) @ np.array(results.params) # params from ols model
quantity_demanded = quantity_demanded[0]
print(quantity_demanded)
我们从方程 8 得到的估计需求量为：
估计需求量函数
现在我们可以构建我们的收入、成本，并将它们结合起来构建我们的利润函数。这里我们生产每单位的成本为 $140 基础，并且每生产一个额外单位的折扣为 $0.0001：
Revenue = price * quantity_demanded
Cost = quantity_demanded * (140 - 0.0001*quantity_demanded) + τ + δ
profit = Revenue - Cost
print(profit)
估计利润函数
将我们的利润绘制为数字广告支出和电视广告支出的函数，*π(*δ、τ)：
现在我们使用 python 按照方程 12 中的公式解决我们的优化问题，利用我们在第一部分和第二部分中学到的优化理论。注意，ρ 的极高值是为了考虑到我们的目标函数值极大，因此我们需要确保惩罚足够大，以避免“跳出”约束。
## Optimization Problem
objective = -profit + λ*(τ + δ - 100_000) - ρ*sm.log((τ-20_000)*(δ-10_000))
symbols = [δ, τ, λ, ρ]
x0 = {δ:20_000, τ:80_000, λ:0, ρ:100000}
results = constrained_newton_method(objective,symbols,x0,iterations=1000)
以及相应的输出：
因此，我们的解决方案是将大约 $66,831 用于数字广告支出，将大约 $33,169 用于电视广告支出。这些数值对应于：
digital_ad = results[δ]
television_ad = results[τ]
quantity = quantity_demanded.evalf(subs={δ:digital_ad,τ:television_ad})
revenue = Revenue.evalf(subs={δ:digital_ad,τ:television_ad})
cost = Cost.evalf(subs={δ:digital_ad,τ:television_ad})
profit = revenue - cost
print(f"Quantity: {int(quantity):,}")
print(f"Total Revenue: ${round(revenue,2):,}")
print(f"Total Cost: ${round(cost,2):,}")
print(f"Profit: ${round(profit,2):,}")
就这样！
结论
在本文中，利润最大化问题绝非完全综合的解决方案。实际上，我们甚至不需要对这么简单的优化问题使用牛顿法！但是，随着优化问题在复杂性和维度上的增加，这在现实世界中非常常见，这些工具变得越来越相关。我们的目标是利用我们在第一部分和第二部分中学到的知识，进行一次有趣的探索，以了解优化理论的无数应用之一。
如果你读到这里，感谢你抽出时间阅读我的文章，并向那些阅读了系列中所有 3 部分的读者表示特别的感谢。我希望你现在对基本的多维优化理论和涉及目标函数约束的扩展有很好的理解。一如既往，我希望你像我写这篇文章时那样享受阅读。如果你对这篇文章和整个系列有任何想法，请告诉我！
附加内容 — 线性回归的数值和解析解
如上所承诺，本节将提供利用牛顿方法解决线性回归问题的代码，并将该结果与解析解（statsmodels 使用的解）进行比较。请回忆，解决线性回归的目标是最小化 残差平方和。也就是说，就矩阵而言，
(A1) 最小化残差平方和
因此，使用我们的牛顿方法函数和框架，我们得到：
# Pull all variables in X and create them as SymPy symbols
variablez = list(df.drop(['quantity_demanded','July'],axis=1).columns)
symbols = []
for i in variablez:
i = sm.symbols(f'{i}')
symbols.append(i)
# Create vectors and matrices of outcome (y), covariates (X), and parameters(β)
y = np.array(df['quantity_demanded'])
X = np.array(df.drop(['quantity_demanded','July'],axis=1))
β = np.array(symbols)
# Specify objective function and starting values
objective = (y - X@β).T @ (y - X@β) # Residual Sum of Squares
β_0 = dict(zip(symbols,[0]*len(symbols))) # Initial guess (0 for all)
β_numerical = newton_method(objective,symbols,β_0)
输出如下：
接下来我们将计算解析解。也就是说，如果我们对等式 (A1*)* 求导并将其设为零并解出 β，我们得到：
(A2) OLS 解析解
编写代码如下 (我们还提供了与 statsmodels 比较的解析标准误差，但不再进一步探讨 — 如果你感兴趣，请参见 OLS 维基百科页面 ）：
# OLS Analytical Solution
β_analytical = np.linalg.inv(X.T @ X) @ X.T @ y
# Compute standard errors
df_residuals = len(X)-len(β_analytical)
σ2 = 1/df_residuals*((y-X@β_analytical).T @ (y-X@β_analytical)) # MSE
Σ = σ2 * np.linalg.inv(X.T @ X)
standard_errors = np.sqrt(np.diag(Σ))
与 statsmodels 的所有结果进行比较：
ols_results = pd.DataFrame()
ols_results['variable'] = variablez
ols_results['β_numerical'] = list(β_numerical.values())
ols_results['β_analytical'] = β_analytical
ols_results['std_err_analytical'] = standard_errors
ols_results['β_statsmodels'] = list(results.params) # from statsmodels code above
ols_results['std_err_statsmodels'] = list(results.bse) # from statsmodels code above
ols_results = ols_results.set_index('variable')
ols_results
附录 — 代码
import sympy as sm
import numpy as np
def get_gradient(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
) -> np.ndarray:
"""
Calculate the gradient of a function at a given point.
Args:
function (sm.core.expr.Expr): The function to calculate the gradient of.
symbols (list[sm.core.symbol.Symbol]): The symbols representing the variables in the function.
x0 (dict[sm.core.symbol.Symbol, float]): The point at which to calculate the gradient.
Returns:
numpy.ndarray: The gradient of the function at the given point.
"""
d1 = {}
gradient = np.array([])
for i in symbols:
d1[i] = sm.diff(function, i, 1).evalf(subs=x0)
gradient = np.append(gradient, d1[i])
return gradient.astype(np.float64)
def get_hessian(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
) -> np.ndarray:
"""
Calculate the Hessian matrix of a function at a given point.
Args:
function (sm.core.expr.Expr): The function for which the Hessian matrix is calculated.
symbols (list[sm.core.symbol.Symbol]): The list of symbols used in the function.
x0 (dict[sm.core.symbol.Symbol, float]): The point at which the Hessian matrix is evaluated.
Returns:
numpy.ndarray: The Hessian matrix of the function at the given point.
"""
d2 = {}
hessian = np.array([])
for i in symbols:
for j in symbols:
d2[f"{i}{j}"] = sm.diff(function, i, j).evalf(subs=x0)
hessian = np.append(hessian, d2[f"{i}{j}"])
hessian = np.array(np.array_split(hessian, len(symbols)))
return hessian.astype(np.float64)
def newton_method(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
iterations: int = 100,
) -> dict[sm.core.symbol.Symbol, float] or None:
"""
Perform Newton's method to find the solution to the optimization problem.
Args:
function (sm.core.expr.Expr): The objective function to be optimized.
symbols (list[sm.core.symbol.Symbol]): The symbols used in the objective function.
x0 (dict[sm.core.symbol.Symbol, float]): The initial values for the symbols.
iterations (int, optional): The maximum number of iterations. Defaults to 100.
Returns:
dict[sm.core.symbol.Symbol, float] or None: The solution to the optimization problem, or None if no solution is found.
"""
x_star = {}
x_star[0] = np.array(list(x0.values()))
# x = [] ## Return x for visual!
print(f"Starting Values: {x_star[0]}")
for i in range(iterations):
# x.append(dict(zip(x0.keys(),x_star[i]))) ## Return x for visual!
gradient = get_gradient(function, symbols, dict(zip(x0.keys(), x_star[i])))
hessian = get_hessian(function, symbols, dict(zip(x0.keys(), x_star[i])))
x_star[i + 1] = x_star[i].T - np.linalg.inv(hessian) @ gradient.T
if np.linalg.norm(x_star[i + 1] - x_star[i]) < 10e-5:
solution = dict(zip(x0.keys(), x_star[i + 1]))
print(f"\nConvergence Achieved ({i+1} iterations): Solution = {solution}")
break
else:
solution = None
print(f"Step {i+1}: {x_star[i+1]}")
return solution
def constrained_newton_method(
function: sm.core.expr.Expr,
symbols: list[sm.core.symbol.Symbol],
x0: dict[sm.core.symbol.Symbol, float],
iterations: int = 100,
) -> dict[sm.core.symbol.Symbol, float] or None:
"""
Performs constrained Newton's method to find the optimal solution of a function subject to constraints.
Parameters:
function (sm.core.expr.Expr): The function to optimize.
symbols (list[sm.core.symbol.Symbol]): The symbols used in the function.
x0 (dict[sm.core.symbol.Symbol, float]): The initial values for the symbols.
iterations (int, optional): The maximum number of iterations. Defaults to 100.
Returns:
dict[sm.core.symbol.Symbol, float] or None: The optimal solution if convergence is achieved, otherwise None.
"""
x_star = {}
x_star[0] = np.array(list(x0.values())[:-1])
optimal_solutions = []
optimal_solutions.append(dict(zip(list(x0.keys())[:-1], x_star[0])))
for step in range(iterations):
# Evaluate function at rho value
if step == 0: # starting rho
rho_sub = list(x0.values())[-1]
rho_sub_values = {list(x0.keys())[-1]: rho_sub}
function_eval = function.evalf(subs=rho_sub_values)
print(f"Step {step} w/ {rho_sub_values}") # Barrier method step
print(f"Starting Values: {x_star[0]}")
# Newton's Method
for i in range(iterations):
gradient = get_gradient(
function_eval, symbols[:-1], dict(zip(list(x0.keys())[:-1], x_star[i]))
)
hessian = get_hessian(
function_eval, symbols[:-1], dict(zip(list(x0.keys())[:-1], x_star[i]))
)
x_star[i + 1] = x_star[i].T - np.linalg.inv(hessian) @ gradient.T
if np.linalg.norm(x_star[i + 1] - x_star[i]) < 10e-5:
solution = dict(zip(list(x0.keys())[:-1], x_star[i + 1]))
print(
f"Convergence Achieved ({i+1} iterations): Solution = {solution}\n"
)
break
# Record optimal solution & previous optimal solution for each barrier method iteration
optimal_solution = x_star[i + 1]
previous_optimal_solution = list(optimal_solutions[step - 1].values())
optimal_solutions.append(dict(zip(list(x0.keys())[:-1], optimal_solution)))
# Check for overall convergence
if np.linalg.norm(optimal_solution - previous_optimal_solution) < 10e-5:
print(
f"\n Overall Convergence Achieved ({step} steps): Solution = {optimal_solutions[step]}\n"
)
overall_solution = optimal_solutions[step]
break
else:
overall_solution = None
# Set new starting point
x_star = {}
x_star[0] = optimal_solution
# Update rho
rho_sub = 0.9 * rho_sub
return overall_solution
通过此 GitHub 仓库访问所有代码： github.com/jakepenzak/Blog-Posts
感谢你阅读我的帖子！我在 Medium 上的帖子旨在探索利用 计量经济学 和 统计/机器学习 技术的实际和理论应用。此外，我还通过理论和模拟提供关于各种方法学的理论基础的帖子。最重要的是，我写作是为了学习！我希望使复杂的话题对所有人稍微更易于理解。如果你喜欢这篇文章，请考虑 关注我在 Medium 上的账号！
优化还是架构：如何破解卡尔曼滤波
为什么即使在神经网络不如卡尔曼滤波器时，它们可能看起来更好——以及如何修复这个问题并改进你的卡尔曼滤波器
·
关注 发表在 Towards Data Science ·6 分钟阅读·2023 年 12 月 26 日
–
本文介绍了我们在 NeurIPS 2023 上发表的 近期论文。代码可在 PyPI 上获取。
背景
卡尔曼滤波器 (KF) 自 1960 年以来一直是一种备受推崇的序列预测和控制方法。尽管在过去几十年中出现了许多新方法，但 KF 的简单设计使其至今仍然是一种实用、稳健且具有竞争力的方法。1960 年的原始论文在过去 5 年内已有12K 次引用。其广泛的应用包括导航、医疗治疗、市场分析、深度学习 甚至登月。
KF 在阿波罗任务中的示意图（图源：作者）
从技术上讲，KF 从一系列噪声观测（例如雷达/相机）中预测系统的状态 x（例如航天器位置）。它估计状态的分布（例如位置估计 + 不确定性）。每个时间步，它根据动态模型 F 预测下一个状态，并根据动态噪声 Q 增加不确定性。每个观测，它根据新的观测 z 及其噪声 R 更新状态及其不确定性。
KF 的单步操作示意图（图源：作者）
卡尔曼滤波器还是神经网络？
KF 的预测模型 F 是线性的，这有些限制。因此，我们在 KF 的基础上构建了一个高端神经网络。这使我们比标准 KF 获得了更好的预测准确性！这对我们来说真是好事！
为了最终确定我们的论文实验，我们进行了几项消融测试。在其中一个测试中，我们完全移除了神经网络，仅仅优化了内部 KF 参数 Q 和 R。当这个优化后的 KF 不仅超越了标准 KF，还超越了我那高端的网络时，你可以想象我脸上的表情！同样的 KF 模型，完全相同的 60 年线性架构，仅通过更改噪声参数的值就变得优越了。
KF、神经 KF (NKF) 和优化 KF (OKF) 的预测误差（图源：我们的论文）
KF 击败神经网络虽然有趣，但对我们的问题来说却是轶事性的。更重要的是方法论的见解：在扩展测试之前，我们差点就宣称网络优于 KF——仅仅因为我们没有正确地比较这两者。
信息 1：为了确保你的神经网络确实比 KF 更好，就像优化网络一样优化 KF。
备注——这是否意味着 KF 比神经网络更好？ 我们确实没有做出这样的普遍性声明。我们的声明是关于方法论的——如果你想比较这两种模型，它们都应当类似地优化。话虽如此，我们确实从经验上展示了在多普勒雷达问题中，尽管存在非线性，KF 可能表现更好。事实上，这让我难以接受，以至于我输了关于我的神经 KF 的赌注，花费了许多周的超参数优化和其他技巧。
优化卡尔曼滤波器
在比较两个架构时，类似地优化它们。这听起来有些微不足道，对吧？实际上，这个缺陷并非我们研究中的独特问题：在非线性滤波文献中，线性 KF（或其扩展 EKF）通常作为比较的基准，但很少被优化。事实上，这背后有一个原因：标准 KF 参数“被认为”已经能够提供最优预测，那么何必再进一步优化呢？
KF 参数Q和 R 的标准封闭形式方程。我们关注的是离线数据同时提供状态{x}和观测{z}的设置，因此协方差可以直接计算。其他确定 Q 和 R 的方法通常用于其他设置，例如没有{x}数据的情况。
不幸的是，封闭形式方程的最优性在实践中并不成立，因为它依赖于一组相当强的假设，而这些假设在现实世界中很少成立。事实上，在简单的经典低维多普勒雷达问题中，我们发现了不少于 4 个假设的违背。在某些情况下，这种违背甚至很难察觉：例如，我们模拟了独立同分布的观测噪声——但使用的是球坐标系。一旦转换到笛卡尔坐标系，噪声就不再是独立同分布的了！
多普勒雷达问题中 4 个 KF 假设的违背：非线性动力学；非线性观测；不准确的初始化分布；以及非独立同分布的观测噪声。（图片由作者提供）
信息 2：不要相信 KF 的假设，因此避免使用封闭形式的协方差估计。相反，按照你的损失函数优化参数——就像任何其他预测模型一样。
换句话说，在现实世界中，噪声协方差估计不再是优化预测误差的代理。这种目标之间的差异会导致意外的异常。在一次实验中，我们用一个oracle KF 替代噪声估计，该oracle知道系统中的准确噪声。这个oracle仍然不如优化后的 KF，因为准确的噪声估计不是期望的目标，而是准确的状态预测。在另一次实验中，当 KF 接收到更多数据时，其表现会恶化，因为它实际上追求的是与均方误差（MSE）不同的目标！
测试误差与训练数据大小的关系。标准 KF 不仅劣于优化后的 KF，而且随着数据的增加而恶化，因为其参数并未设置以优化期望的目标。（图像摘自我们的论文）
那么，如何优化 KF 呢？
在标准噪声估计方法用于 KF 调整的背后，是将 KF 参数视为噪声的代表。这种观点在某些情况下是有益的。然而，如上所述，为了优化，我们应该“忘记”KF 参数的这一角色，只将它们视为模型参数，其目标是损失最小化。这种替代观点也告诉我们如何进行优化：就像任何序列预测模型，例如 RNN！给定数据，我们只需进行预测，计算损失，反向传播梯度，更新模型，并重复进行。
与 RNN 的主要区别在于，参数Q和R以协方差矩阵的形式出现，因此它们应该保持对称且正定。为此，我们使用 Cholesky 分解将Q写成LL**的形式，并优化L的条目。这确保了Q在优化更新后仍然保持正定。这个技巧同时适用于Q和R*。
OKF 训练过程的伪代码（摘自我们的论文）
在我们所有的实验中，这种优化过程被发现既快速又稳定，因为参数数量比典型神经网络少几个数量级。虽然训练过程容易自己实现，但你也可以使用我们的 PyPI 包，如这个示例所示 😃
总结
如下图所示，我们的主要信息是 KF 假设不可靠，因此我们应该直接优化 KF——无论我们是将其作为主要预测模型，还是作为与新方法比较的参考。
我们的简单训练过程可以在 PyPI 上找到。更重要的是，由于我们的架构与原始 KF 保持一致，任何使用 KF（或扩展 KF）的系统都可以通过重新学习参数轻松升级为 OKF——无需在推理时增加任何代码行。
我们的范围和贡献总结。由于 KF 假设经常被违反，因此必须直接优化 KF。 （图片来自我们的论文）
优化数据仓库存储：视图与表
原文：
towardsdatascience.com/optimize-data-warehouse-storage-with-views-and-tables-659dd588345d
表与视图的区别及其使用方法
·发布于 Towards Data Science ·阅读时间 6 分钟·2023 年 3 月 24 日
–
照片由 Sophia Baboolal 提供，来源于 Unsplash
随着现代数据技术堆栈的兴起，许多公司正在将其数据库从本地迁移到云端。他们开始利用像 Snowflake、Redshift 和 DuckDB 这样的数据仓库工具，以充分发挥云的所有优势。
尽管这些数据仓库通常帮助较小的公司节省开支，但云上的计算成本很容易累积。因此，优化存储和计算成本至关重要。这意味着你需要了解最佳的数据存储方式，以便数据团队能够以具有成本效益的方式使用这些数据。
在本文中，我们将讨论视图与表的区别、数据仓库中存在的不同类型视图以及每种视图的使用场景。阅读完本文后，你应该能够确定存储不同数据集的最佳选项，同时节省成本。
什么是视图？
视图是一个定义的查询，位于表的顶部。与表不同，它不存储实际数据。它始终包含最新的数据，因为每次查询时都会重新运行。而表的内容仅与上次创建或更新时一样，无论你什么时候查询。
视图主要有两种类型：非物化视图和物化视图。
非物化视图
非物化视图是人们通常想到的视图。这种类型的视图仅在实际查询时运行，否则不会存储在数据库中。
非物化视图很棒，因为它们不占用存储空间，这意味着你不必担心为大量存储付费。它们也只有在需要时才运行，从而节省计算资源费用。这意味着，如果源表几个月或几周都不需要，你无需为其维护付费。只有当分析师或分析工程师重新开始使用该表时，你才需为其付费。
最棒的部分？非物化视图仍然具有与表相同的所有功能！如果需要，你可以在它们上面执行连接、聚合和窗口函数。
不幸的是，就像其他所有事物一样，所有优点背后总会有缺点。非物化视图不适合处理大量数据和复杂逻辑，因为每次查询视图时，这些逻辑都会被执行。
例如，我通常将所有的源数据表创建为非物化视图，这些视图引用我的原始数据。这些是简单的 SELECT 语句，包含基本函数，如列重命名、数据类型转换和数据清理。由于它们的底层逻辑简单，所以每当我查询这些源表时，它们运行得很快。
如果我创建包含连接和窗口函数的复杂数据模型作为视图，可能在查询时这些视图永远不会加载。或者它们需要极其长的时间！显然，这不是理想的情况。你将会花费更多的计算能力来运行视图查询，而不是将视图创建为表。
记住：非物化视图非常适合使用，但仅当创建它们的逻辑是简单的 SELECT 语句时。
物化视图
物化视图在我们讨论的两种视图中较为少见。物化视图的行为更像表。它们查询速度更快，被认为比非物化视图更易于访问。并且，就像表一样，它们在数据仓库中占用更多的存储空间，需要更多的计算资源。这意味着它们是两种视图中更昂贵的选项。
你不会经常需要使用它们。事实上，我从未遇到过一个合理的使用场景。根据Snowflake 的文档，你应该仅在以下所有条件都满足时使用物化视图：
-
视图的结果被频繁使用
-
驱动视图的查询使用了大量资源
-
视图变化频繁
这三者同时适用于你的基础/暂存、中间和核心 dbt 模型的情况非常少见。基础/暂存模型不会消耗大量资源，而中间和核心数据模型变化不频繁。当然，总会有例外情况，但我还没遇到过这种情况。
这些视图在数据建模中的使用方式
如果你是分析工程师，你可能会想知道未物化视图和物化视图在数据建模中如何使用。让我们深入了解一下 dbt 基础（或阶段）模型以及核心模型。
dbt 基础模型
dbt 基础模型作为视图存在于你的原始数据之上。它们被创建为未物化视图，以保持原始数据的完整性，同时利用适当的命名约定和公司标准。这些模型中的代码是从原始数据中直接读取的基本 SQL 选择语句，通过像 Airbyte 这样的摄取工具进行 ELT。一个典型的基础模型如下：
select
ad_id AS facebook_ad_id,
account_id,
ad_name AS ad_name_1,
adset_name,
month(date) AS month_created_at,
date::timestamp_ntz AS created_at,
spend
from {{ source('facebook', 'basic_ad')}}
如果你查看 dbt 中此文件的底层逻辑，它实际上会在 Snowflake（我首选的数据仓库）中编译为如下所示：
create or replace view data_mart_dev.base.base_facebook_ads
as (
select
ad_id AS facebook_ad_id,
account_id,
ad_name AS ad_name_1,
adset_name,
month(date) AS month_created_at,
date::timestamp_ntz AS created_at,
spend
from raw.facebook.basic_ad
);
由于你仅使用基本的日期函数和重命名列，视图在按需查询时仍然很快。这反过来节省了你原本会用来保存几乎相同的原始数据副本的存储空间。
dbt 核心模型
在 dbt 中，你的核心模型比基础模型更复杂，通常包含多个 CTE、联接和窗口函数。虽然你可能有特定的用例来将这些创建为物化视图，但你很可能会将这些创建为数据仓库中的表。表对于处理复杂的转换非常理想，如果将其存储为视图，会需要较长的运行时间。
这是我一个核心数据模型的代码示例：
with
fb_spend_unioned AS (
select created_at, spend, 'company_1' AS source from {{ ref('base_fb_ads_company1')}}
UNION ALL
select created_at, spend, 'company_2' AS source from {{ ref('base_fb_ads_company2')}}
),
fb_spend_summed AS (
select
month(created_at) AS spend_month,
year(created_at) AS spend_year,
created_at AS spend_date,
sum(spend) AS spend
from fb_spend_unioned
where spend != 0
group by
created_at,
month(created_at),
year(created_at)
)
select * from fb_spend_summed
当在 Snowflake 中编译为 SQL 时，代码将如下所示：
create or replace table data_mart_dev.core.fb_spend_summed
as (
with
fb_spend_unioned AS (
select created_at, spend, 'company_1' AS source from {{ ref('base_fb_ads_company1')}}
UNION ALL
select created_at, spend, 'company_2' AS source from {{ ref('base_fb_ads_company2')}}
),
fb_spend_summed AS (
select
month(created_at) AS spend_month,
year(created_at) AS spend_year,
created_at AS spend_date,
sum(spend) AS spend
from fb_spend_unioned
where spend != 0
group by
created_at,
month(created_at),
year(created_at)
)
select * from fb_spend_summed
);
请注意，这里创建的是一个 Snowflake 中的表而不是视图。这对于任何将直接用于 BI 工具的数据是理想的，大多数核心数据模型都是这样。它们可以按需轻松查询，无需运行底层逻辑。这确保了快速的仪表板，能够让利益相关者信赖。
结论
视图和表在数据仓库中存在不同的原因。视图不存储实际数据，可以作为一个工具来节省开支，通过简单查询在其他表之上运行。表应当用于存储由更复杂逻辑生成的数据，确保性能和可用性始终较高。
正确使用时，非物化视图是在 Snowflake 中节省开支的好工具，同时不牺牲性能。我强烈推荐在 dbt 中将其用于你的基础模型，以创建符合你设定的公司标准的高质量数据。而且，不要忘记在核心 dbt 模型中使用表。性能提升是值得更高成本的！
欲了解更多分析工程的信息，订阅我的免费每周通讯，我会分享学习资源、教程、最佳实践等。
查看我的第一本电子书，分析工程的基础知识，这是一本关于如何开始从事分析工程角色的全方位指南。
优化浏览分类法
原文：
towardsdatascience.com/optimizing-browsing-a-taxonomy-30596bcdbd9d
问题、模型和实施问题
·发表在Towards Data Science ·阅读时间 20 分钟·2023 年 1 月 4 日
–
图片来自mcmurryjulie，pixabay
这里的主要用例是jagota-arun.medium.com/interactive-and-adaptive-menu-ordering-agent-bb447c58b3af（不是必读材料。）
我们专注于建模。它是什么？它如何与推荐系统相关？然后是实施。就像面向对象设计一样。类等。
面向对象的阐述风格将序列化、模块化、可组合性和可扩展性的好处扩展到写作中。这里非常合适。
本文最重要的概念是分类法。一切都发生在它上面。用户浏览、用户行为（在叶子上）、建模和评分，以动态重新排列分类法，优化用户的浏览体验。
那么让我们从定义它开始。分类法是一个选择树，具体的感兴趣项目位于其叶子上。我们最终假设用户关心的是到达感兴趣的叶子。内部节点主要用于导航。我们通过用户是否在叶子上采取行动来学习。
一个例子是餐厅菜单。菜单项位于叶子上。内部节点是菜单上的各种类别：饮料、开胃菜、主菜等。用户的一个动作是点选特定的菜单项。如果用户想点饮料，她会进入饮料子树。
在视频观看场景中，用户的行为可能是观看视频。
我们将用户与分类法的互动建模为一个会话，其中用户浏览树的某些部分。当用户到达一个合适的叶子时，她可以选择行动或不行动。
我们的主要兴趣是从用户在特定叶子上的行为中学习偏好，以便在树中的某些节点上表现出偏好。我们希望在叶子节点和内部节点上都能学习到偏好。
这与常规推荐系统有何关系？
通常，我们期望推荐系统主要推荐新的项目。我们在这里的重点是简化用户对感兴趣项目的导航。这可能包括用户之前已操作过的项目。例如餐厅设置。
另一个需要注意的是，推荐系统通常基于用户-物品矩阵操作，而不是基于分类法。
面向对象的阐述
我们最重要的类将是 UserInteractionsWithTaxonomy。其角色应从其名称和我们之前的讨论中明确。
这个类会被实例化以建模特定用户与特定分类法的互动。实例化时需要提供分类法。调用者负责将实例化的模型与适当的用户关联。如下面所示。这个类会为分类法中的每个节点维护一个计数器。这些计数器在新实例中会被初始化为零。
john_doe_interactions_with_specific_restaurant =
UserInteractionsWithTaxonomy(taxonomy_of_restaurant)
“用户”可以是模型创建者想要的任何人。一个人或一整个群体。
关键方法
UserInteractionsWithTaxonomy 将支持一个方法 act。这个方法将建模用户在特定叶子上的一次行为实例。例如，点餐时选择某个特定的菜品。
在我们版本的 act 方法中，每次调用都会递增叶子及其所有祖先的计数。为了高效地做到这一点，这个类可能会维护分类法中所有节点的父节点弧。
请注意，我们选择仅从 act 行为中学习用户偏好。我们可以将此推广到包括用户点击树的内部节点的反馈。但在这里我们不会追求这一点，因为我们认为 act 行为比内部节点点击更能预测用户偏好。读者可以自行尝试。
第二个方法将模拟用户在树上的某个节点上点击。这将让我们追踪用户当前在树上的位置，以便我们可以在下面展示适当的选择。
我们将在一个方法 choices() 中表达这些选择，该方法将返回当前用户在树上所在节点的子节点的得分。更精确地说，它将返回一个 [child, score] 对列表，其中 child 是 self.current_node 的子节点，score 是其得分。这个列表按得分的非递增顺序排列。
这是一个示例流程。
john_doe_interactions_with_specific_restaurant.click("Drinks")
choices = john_doe_interactions_with_specific_restaurant.choices()
# User picks 'Coca Cola' from choices and acts on it.
# ('act' is 'order here'.)
john_doe_interactions_with_specific_restaurant.act('Coca Cola')
让我们一起来看看这个类的内部细节。
class UserInteractionsWithTaxonomy:
def _init_(self, taxonomy):
self.taxonomy = taxonomy
self.current_node = None
def click(self, node):
self.current_node = node
def act(self, leaf, pc = 1):
# increment the leaf's counter by pc
# increment the counters of all of the leaf's ancestors, each by pc
def choices(self):
res = {}
for each child of self.current_node:
res[child] = self.score(child, self.current_node)
return res
def score(self,child,parent):
return self.counter[child]/self.counter[parent]
请注意 act 方法中的参数 pc。我们稍后会解释它。目前，假设 pc 等于 1，即其默认值。因此，“递增 … by pc” 意味着 “递增 … by 1” 或简称 “递增 …”。
方法score(.)返回P(child|parent)，即用户在parent位置时访问child的概率。这只是子项的计数除以父项的计数。
冷启动问题
假设用户访问了一家新餐馆。由于用户在该餐馆没有订单历史记录，因此不能使用评分方法。因此，我们必须等用户点几个项目后才能评分。这个场景被称为用户冷启动问题。
类似地，还有一个分类冷启动问题。在餐馆的场景中，这对应于一个全新的餐馆开张。尚未与任何用户互动。
实际上，我们当前评分中隐藏了更细致的冷启动问题。我们的评分函数是
self.counter[child]/self.counter[parent]
想象一下一个餐馆，用户有订单历史记录，但在特定类别下没有，例如甜点。父类别 = 甜点的计数将为 0，因此其子类别的得分也无法计算。
现在想象一下，用户在某个类别中下单，但未在特定子类别中下单。例如，在饮料类别中点了一种特定的饮料。子项的得分将为 0。虽然 0 是一个有效的得分，但可能不是我们所希望的。
总结来说，未曾被点过的项目将无法评分或得分为 0。
冷启动问题在推荐系统的环境中经常出现，见[4]。在这种环境下，分类冷启动问题对应于项目冷启动问题。
缓解冷启动问题
缓解所有冷启动问题的一个简单方法是将所有节点上的计数器初始化为一个称为pc的正数。虽然它需要是正数，但可以小于 1。
在我们的设置中，我们将支持一个类似但更强大的机制。无论是从建模还是编码的角度来看，都不会更复杂。
我们将为类UserInteractionsWithTaxonomy添加一个方法act_randomly(n, pc)。
class UserInteractionsWithTaxonomy:
…
def act_randomly(self, n. pc):
for i in range(n):
leaf = leaves.sample_randomly()
self.act(leaf, pc)
现在回到用户第一次访问餐馆的场景。假设我们调用act_randomly(n = 100, pc = 0.01)。这样做的效果就像用户在这家餐馆点了 100 次，每次随机选择菜单上的一项。为了考虑pc，我们应该将这些订单称为“分数订单”。
为什么要使用pc？也就是说，为什么不设置为 1？因为在某些环境中，我们可能希望单个实际订单比单个虚拟订单更重要。在这种情况下，我们会将pc设置为远小于 1 的值。在其他环境中，我们可能希望相反的行为。只有当用户订单达到一定的最小次数时，我们才会认为它是重要的。在这种情况下，我们会将pc设置为大于 1 的值。
act_randomly方法可以缓解所有上述冷启动问题。不过，我们通常可以在用户冷启动问题上做得更好。
使用分类模型加热用户冷启动
考虑一个有很多用户互动历史的餐厅。这些历史互动可以用来为该餐厅形成一个模型。这个模型，我们称之为餐厅模型，实际上是UserInteractionsWithTaxonomy的一个实例，其中用户是一个集体而非个人。
为了在特定用户的交互中利用分类模型，我们将重构UserInteractionsWithTaxonomy类如下。
class UserInteractionsWithTaxonomy:
def __init__(self, taxonomy, taxonomy_model = None):
…
self.taxonomy_model = taxonomy_model
…
def act_on_this(self, leaf, pc = 1):
# Previous version of act(.) renamed
def act(self, leaf):
self.act_on_this(leaf)
if self.taxonomy_model is not None:
self.taxonomy_model.act_on_this(leaf)
def act_randomly(self, n. pc):
for i in range(n):
…
self.act_on_this(leaf, pc)
def score(self,child,parent, pcm = 1):
if self.taxonomy_model is None:
return self.counter[child]/self.counter[parent]
else:
s = self.taxonomy_model.score(child, parent)
pc = pcm*s
child_adjusted_count = self.counter[child] + pc
parent_adjusted_count = self.counter[parent] + pcm
return child_adjusted_count/parent_adjusted_count
让我们解释一下。在实例化UserInteractionsWithTaxonomy时，我们提供另一个UserInteractionsWithTaxonomy的实例，它将用作分类模型，如下所示。
restaurant_model = UserInteractionsWithTaxonomy(taxonomy_of_restaurant)
user_restaurant_models = {}
for user in users_of_this_restaurant:
user_restaurant_models[user] = UserInteractionsWithTaxonomy(
taxonomy_of_restaurant, restaurant_model)
在上面的示例中，所有用户餐厅模型共享相同的（通用）餐厅模型。
对任何用户餐厅模型上的act，即使是未来的模型，也会对共享餐厅模型执行相同的操作。
为了干净地实现这一点，我们将之前的act版本重命名为act_on_this，以便可以在其上构建新的act。由于这一名称更改，我们还不得不对act_randomly(.)进行更改。
然后我们重构了方法score(.)，使得如果实例附加了分类模型，它也会在评分中使用该模型。
在score(.)的重构版本中发生了很多事情，所以我们来详细说明一下。考虑在用户特定模型中对某个孩子与某个父母进行评分。首先，我们调用
self.taxonomy_model.score(child, parent)
这给我们提供了P(child|parent, generic taxonomy model)。接下来，我们假设对父母进行pcm次访问，其中pcmP(child|parent, generic taxonomy model)预计会以child*结束。因此，我们将这些想象中的访问次数加到用户的实际访问次数上。
在上面的段落中，术语“访问节点”仅用于说明。在实际情况中，对一个节点的访问对应于在该节点下的叶子上的act实例。
让我们看一个数字示例。假设pcm是 10。假设P(child|parent, generic taxonomy model)等于⅕。对于这个孩子，pc将是 2。如果特定用户从未访问过父母，P(child|parent, user model for this taxonomy)将等于⅕，即P(child|parent, generic taxonomy model)的值。随着用户开始访问parent，用户随后访问的孩子将开始影响P(child|parent, user model for this taxonomy)，如果用户的行为偏离，可能会将其从P(child|parent, generic taxonomy model)中移开。
建模用户偏好
考虑这种情况。一个全新的餐厅刚刚开张。由于没有历史记录，因此没有餐厅模型。
一个用户走进来。当然，我们可以使用act_randomly从菜单中随机采样，如前所述。我们能做得更好吗？
我们应该将餐馆模型实例与每一个访问餐馆的新用户关联起来。这样，至少餐馆模型会比任何一个个人学习得更快。任何突出显示的菜单项，即在用户中被频繁点餐的项，会更快地被呈现出来。我们还能做得更好吗？
假设用户过去访问过其他餐馆。如果我们能够对这些访问建模并学习用户的偏好，我们就可以将这些偏好应用到用户首次访问这家新餐馆时。例如，某些词汇如curry的偏好，或某些价格范围的偏好。
虽然在前一段中我们使用了餐馆冷启动问题来激励建模用户偏好，但解决这个问题涉及到复杂的知识转移问题。这是因为用户之前访问的餐馆的分类法可能与当前餐馆的分类法不同。所以我们将推迟讨论这个用例到另外一篇文章中。
也就是说，将用户偏好纳入我们的模型可以帮助未来用户与相同分类法的互动。
比如，当有很多菜品时。用户可能会错过一些好的菜品。一些餐馆的菜单上有超过四十种菜品。在其他情况下，分类法甚至可能更大。例如，Netflix、Amazon Prime 或其他各种流媒体视频服务上的视频。
或者，当相同的词出现在多个项目中时。这甚至适用于小型分类法。想象一下访问一家泰国餐馆，该餐馆有几个主菜中包含了词汇 curry。假设用户点了其中一道菜。下次，用户会看到其他咖喱菜品被显著展示。这很有道理，对吧？
词袋模型
既然我们决定对用户偏好进行建模，就让我们开始吧。我们将把它们建模为词袋。
词袋模型是一个词频对的集合。一个例子是 { beer :10, wine :2 }。我们可以将这个例子解释为用户偏好啤酒而非葡萄酒。
虽然我们在这里使用了“词”这个术语，但我们可以根据需要定义它。例如，我们可以把术语“wheat beer”表示为“wheat-bear”。
我们将把BagOfWords建模为一个类。我们将用一个哈希来表示词袋，即 Python 的Dict。键是单词，值是它们的频率。
我们将如何具体使用BagOfWords？在实例化UserInteractionsWithTaxonomy时，我们还将创建一个BagOfWords实例来建模用户的关键词偏好。当用户对某个叶子进行操作时，我们将把该叶子名称的词袋添加到代表用户偏好的词袋中。为此，BagOfWords 类需要支持一个“add”方法，用于将提供的词袋添加到实例的词袋中。
这是我们的初始版本的 BagOfWords。
class BagOfWords:
def _init_(self):
# Create self.bow
def add(self, bow):
for word in bow:
self.bow[word] += bow[word]
这是用户与分类法交互的版本，重构以建模她在分类法中的交互中的用户偏好。
class UserInteractionsWithTaxonomy():
def _init_(self, taxonomy):
…
self.user_likes = BagOfWords()
def act(self, leaf):
…
self.user_likes.add(leaf.name.bow)
使用用户偏好
现在我们完成了首个版本的用户偏好建模和训练，让我们讨论一下如何使用它。
我们从回顾UserInteractionsWithTaxonomy中的score方法开始。它对特定子项进行评分，以区分给定父项的各种子项，用于导航目的。
现在假设我们想要增强评分功能，以便除了订单数量外，还使用用户偏好。
在我们当前的设计中，用户偏好被建模为在用户与分类法交互实例中维护的单一词袋。鉴于我们打算如何使用评分函数，即区分一个父项的各种子项，如果每个父项拥有自己的用户偏好会更自然。这意味着分类法中的每个内部节点都应该有自己的用户偏好。
这将使我们可以将 UserInteractionsWithTaxonomy.score 重构如下。
class UserInteractionsWithTaxonomy:
…
def score(self,child,parent, pcm = 1):
…
if it makes sense to score off user likes and not order frequencies:
return self.user_likes[parent].score(self.user_likes[child])
…
我们不会详细说明“是否有意义基于用户喜欢而非订单频率进行评分”这一点。这需要一些讨论。
从单一用户偏好模型转向节点特定的用户偏好模型乍看之下可能令人望而生畏。经过进一步考虑，首先，这种变化并不像最初看起来那样令人生畏；其次，它实际上可以使评分更加准确。
详细说明一下，实际中许多分类法会很小。例如，即使是一个全方位服务的餐厅的分类法，通常也只有四个内部节点，不包括根节点：饮料、开胃菜、主菜和甜点。
详细说明一下，某些词可能跨越节点边界。在某个餐厅中，例如，三文鱼这个词出现在开胃菜的名称中，也出现在主菜的名称中。假设从某个用户的订单历史中，我们发现她喜欢三文鱼开胃菜但不喜欢三文鱼主菜。当她寻找开胃菜时，我们应该显著展示三文鱼开胃菜。当她寻找主菜时，我们则不应显著展示三文鱼主菜。
既然我们决定继续进行这种建模增强，我们需要重构*init()以构造节点特定的BagOfWords实例，并适当地更新act*(.)。以下是重构后的版本。
class UserInteractionsWithTaxonomy:
def _init_(self, taxonomy):
…
self.user_likes = {}
for each internal node in the taxonomy:
self.user_likes[node] = BagOfWords()
def act(self, leaf):
…
for each node that is an ancestor of leaf:
self.user_likes[node].add(leaf.name.bow)
何时基于用户偏好而非订单频率进行评分
现在让我们讨论一下
if it makes sense to score off user likes and not order frequencies:
return self.user_likes[parent].score(self.user_likes[child])
部分。
当自我.计数器[父项]的值足够大时，基于订单频率来进行评分比基于用户偏好更好。这是因为实际订单通常应该比用户偏好具有更高的权重。某个特定项目中可能有些东西未被其名称捕捉到，从而解释了用户为何频繁订购它。或者避免它。也许是口味。
当self.counter[parent]的值不够大时，我们有几个选择。我们可以基于分类模型和用户的订单历史来打分。或者我们可以根据用户的偏好来打分。
用简单的语言来说，上一段的选项可以理解为以下内容。
对某个项目或类别与其父类别进行评分，依据
-
其他人在这家餐厅的体验。
-
用户自己在这个父类下的体验，可能结合其他用户的体验，点了与这个词有一些相同的词的项目。
让我们用一个例子来说明这两者。假设一个用户第一次访问某个泰国餐馆。如果用户知道某道主菜非常受欢迎，她可能会有兴趣点这道菜。如果用户知道餐馆中的主菜中含有咖喱这个词很受欢迎，那么所有这些主菜都应该为这个用户打高分。
BagOfWords 中的评分方法
首先，让我们提醒自己我们在UserInteractionsWithTaxonomy.score(.)中添加的相关部分：
if it makes sense to score off user likes and not order frequencies:
return self.user_likes[parent].score(self.user_likes[child])
我们看到我们需要向BagOfWords添加一个score()方法。调用BagOfWords.score(bag_of_words)将把被调用的对象解释为父对象，把传递给它的词袋解释为某个子对象。
我们之所以提到这一点，是因为这里存在固有的不对称性。我们不是试图以对称的方式对两个词袋进行评分。相反，我们寻求评分一个特定子词袋在其父词袋中的适配程度。
让我们看一个例子来解释这一点，同时给出一些对我们评分函数的具体要求。
考虑一个有啤酒类别并且有多个项目的餐厅。假设一个用户访问过餐厅几次，每次都点了淡啤酒。假设淡和啤酒这两个词在这些点过的项目中明确出现。用户在啤酒类别下的词袋现在将包含这两个词。
假设在啤酒菜单中有一个名为Coors Light的项目，但用户在之前的访问中没有点过这个项目。这个项目应该得高分，因为它是一种淡啤酒。即使用户还没有点过包含Coors这个词的啤酒。
现在考虑啤酒类别下的第二个项目，名为Coors Beer。它的评分不应该像Coors Light那么高。为什么呢？因为啤酒这个词在啤酒类别中的信息量不如淡。假设这两个词在用户的啤酒类别词袋中的频率相同，那么淡应该胜出。
现在，让我们朝着一个能满足我们需求的评分函数努力。
我们的第一个决定是仅对那些在子词袋和父词袋中都存在的词进行评分。对于我们的两个例子，如下所示。
{ coors, light } and { light, beer } ⇒ light
{ coors, beer } and { light, beer} ⇒ beer
我们的第二个决定是独立并且加法地评分这个交集中的词。
我们的第三个决定是对交集中的任何单词 w 进行如下评分
f(w, child)*log P(w|parent, user)/P(w|parent, uniform user)
这里 f(w, child, user) 是用户的子词袋中 w 的频率。
P(w|parent, user) 是 f(w, parent, user)/sum_{w’ in parent} f(w’, parent, user)
P(w|parent, uniform user) 类似，只是用户是我们所说的“均匀用户”。均匀用户被假定为对 child 下的任何项都有相等的订购概率。
在关联规则挖掘中，P(w|parent, user)/P(w|parent, uniform user) 是所谓的提升 [3] 的一种版本。
让我们看看这些在我们的示例中是如何工作的。下图展示了这个用户对单词 light 和 beer 的评分。为了简化表达，省略了术语 parent。
log P(light|user)/P(light|uniform user)
log P(beer|user)/P(beer|uniform user)
P(light|user) 和 P(beer|user) 各为 ½。假设菜单上包含单词 beer 的啤酒比包含单词 light 的啤酒多。P(light|uniform user) 可能会显著小于 P(beer|uniform user)。因此，log P(light|user)/P(light|uniform user) 可能会显著大于 log P(beer|user)/P(beer|uniform user)。这是预期中的结果。
让我们在伪代码中查看完整的得分函数。实现中的一个重要细节是评分将分布在两个类之间：BagOfWords 和 UserInteractionsWithTaxonomy。首先，让我们查看两个类中的实际更改，然后再解释原因。
让我们从重构后的 BagOfWords 开始。由于更改并不是局部的，我们将在这里全面展示类的新版本。
class BagOfWords:
def _init_(self):
# Create self.bow
n = 0
def add(self, bow):
for word in bow:
self.bow[word] += bow[word]
n += bow[word]
def p(self, w):
return float(self.bow[w])/n
def score(self, user_bow):
words_in_common = set(user_bow.keys()) and set(self.bow.keys())
score = 0
for w in words_in_common:
score += user_bow[w]*log(self.p(w))
return score
现在让我们引入重构后的 UserInteractionsWithTaxonomy 类。
class UserInteractionsWithTaxonomy:
def __init__(self, taxonomy, taxonomy_model = None, uniform_acts = None):
…
self.taxonomy_model = taxonomy_model
self.uniform_acts = uniform_acts
…
def score_using_user_preferences(self, child, parent):
user_score = self.user_likes[parent].score(self.user_likes[child])
uniform_score = self.uniform_acts.user_likes[parent].
score(self.user_likes[child])
return user_score - uniform_score
def score(self,child,parent, pcm = 1):
…
if it makes sense to score off user likes and not order frequencies:
return self.score_using_user_preferences(child, parent)
…
扩展用户喜好模型
到目前为止，我们已经将用户喜好建模为词袋。在许多用例中，分类法的叶子上可能还有其他属性，这些属性可能会影响用户的喜好。例如 price。
为了超越词袋的模型，将用户喜好建模为一个独立的类 UserLikes 是合理的。我们还需要重构 UserInteractionsWithTaxonomy 以使用 UserLikes。
首先，让我们回顾一下 BagOfWords 在 UserInteractionsWithTaxonomy 中的使用情况，因为这是我们需要推广的内容。
class UserInteractionsWithTaxonomy:
def _init_(self, taxonomy):
…
…
self.user_likes[node] = BagOfWords()
def act(self, leaf):
…
for each node that is an ancestor of leaf:
self.user_likes[node].add(leaf.name.bow)
def score_using_user_preferences(self, child, parent):
user_score = self.user_likes[parent].score(self.user_likes[child])
uniform_score = self.uniform_acts.user_likes[parent].
score(self.user_likes[child])
return user_score - uniform_score
很明显，我们需要进行以下更改：
self.user_likes[node] = BagOfWords() ⇒
self.user_likes[node] = UserLikes()
self.user_likes[node].add(leaf.name.bow) ⇒
self.user_likes[node].add(leaf)
假设这些更改已经完成。
现在让我们讨论一下对 UserLikes 类的影响。
它需要支持一个方法 add(leaf)，该方法从 leaf 中提取所需的所有属性，并将它们添加到当前的用户喜好模型中。
它需要支持一个方法 score，以根据 self（即父节点）对任何给定子节点的用户喜好进行评分。
我们几乎准备好查看UserLikes类了。在此之前，让我们建模一个额外的属性，因为这就是我们首先构建此类的原因。
让我们从价格开始。首先，我们需要讨论其建模方面。
价格基础建模
想象一下访问一家餐馆，其菜单如下所示。这里只显示了部分。对于这些项目，仅显示价格。
饮料：$2, $3, $2, $8, $9, … 主菜：$8, $9, $22, $28, $9, $8, …
“…”是为了表示每个类别中还有许多其他项目。
加粗的项目是特定用户在这家餐馆过去订购的项目。由此，我们可以推断出至少在这家餐馆中，这位用户偏好每个类别中相对较贵的项目。（当然，可能还有其他因素，但这就是我们从这些数据中可以得出的结论。）
这些信息可以用来在用户下次访问这家餐馆时更显著地展示某些项目。也许以不同的顺序展示，如下所示。
饮料：$8, $9, $2, $3, 2 , … ∗ 主菜 ∗ ： 2, … *主菜*： 2,…∗主菜∗：22, $28, $8, $9, $9, $8, …
价格详细建模
我们可以像处理项目名称中的单词一样处理价格。这很好地融入了我们当前的设计。每个内部节点都可以有自己的价格包，就像它有自己的单词包一样。
节点特定的价格包比节点特定的单词包更有动机。这是因为价格在预测的类别上比单词更不具体。用户可能不会在饮料上花费$8，但可能会在主菜上花费这么多。类别特定的价格包将允许我们建模这种区分。
现在，我们将价格包建模为BagOfWords的一个实例。这将满足我们基本的需求。即，对于每个类别，我们记录用户在该类别中订购的项目价格及其频率。
UserLikes
现在是伪代码部分。
class UserLikes:
def __init__(self):
self.bow = BagOfWords()
# Additional attributes
self.prices = BagOfWords()
def add(self, leaf):
self.bow.add(leaf.name.bow)
self.prices.add(leaf.price.as_bow())
def score(self, user_likes):
# This needs some modeling discussion first.
注意leaf.price.as_bow。 “as_bow”是指价格（一个标量）已被表示为一个包含单一键（价格）及其频率为 1 的单词包。
UserLikes.score(.)
考虑一下这个。
class UserLikes:
…
def score(self, user_likes):
return self.bow.score(user_likes.bow) +
self.prices.score(user_likes.prices)
注意到自我是父级的UserLikes，在score(.)的参数中提供的user_likes是特定子级的，通过观察两次调用score(.)中的发生情况，我们可以将其表示为
sum_{w appears in both parent.user_likes.bow and in child.user_likes.bow}
child.user_likes.bow[w]*log(parent.user_likes.bow.p(w)
+
sum_{p appears in both parent.user_likes.prices and in child.user_likes.prices}
child.user_likes.prices[p]*log(parent.user_likes.prices.p(p)
现在考虑到
UserInteractionsWithTaxonomy.score_using_user_preferences(child, parent)
它调用两次，一次如上所述，一次对不同的用户进行均匀交互。
因此，整体评分函数可以表示为
(
sum_{w appears in both parent.user_likes.bow and in
child.user_likes.bow}
child.user_likes.bow[w]*log(parent.user_likes.bow.p(w)
+
sum_{p appears in both parent.user_likes.prices and in
child.user_likes.prices}
child.user_likes.prices[p]*log(parent.user_likes.prices.p(p)
)
-
(
sum_{w appears in both parent.uniform_user.user_likes.bow and in
child.uniform_user.user_likes.bow}
child.uniform_user.user_likes.bow[w]*log(parent.uniform_user.user_likes.bow.p(w)
+
sum_{p appears in both parent.uniform_user.user_likes.prices and in
child.uniform_user.user_likes.prices}
child.uniform_user.user_likes.prices[p]*log(parent.uniform_user_likes.prices.p(p)
)
让我们在一个稍作改进的场景中尝试理解这一点。假设有一家餐厅的分类为啤酒。一个特定的用户已经两次光顾这家餐厅，每次点的都是一些淡色啤酒。（两次点的啤酒不同。）假设“淡色”和“啤酒”这两个词显式出现在所点啤酒的名称中。
现在进行改进。假设两款点的啤酒分别定价为$5 和$6。现在考虑菜单上另外两款啤酒。两者均不包括在用户之前点过的两款啤酒中。
Coors Light $5
? Light $8
上述评分函数会给Coors Light更高的分数。这是因为虽然两者都包含“light”这个词，但第二款的价格显著高于用户之前点的啤酒。
我们可以通过替换进一步改进这一点
p appears in both parent.user_likes.prices and in child.user_likes.prices
由
p is in the price ranges of both parent.user_likes.prices and of child.user_likes.prices
我们不会在这篇文章中再多说了。
总结
在这篇文章中，我们深入探讨了优化用户在分类体系中的浏览，以便在叶子节点中找到感兴趣的项目。我们介绍了一种从用户操作的项目反馈中学习的算法。（在餐厅环境中，即用户实际点的项目。）
反馈被推送到分类体系的上层，以便内部节点可以从中受益。这使得这些节点可以被更（或更少）优先显示，以便用户在未来更容易找到感兴趣的项目。
“感兴趣的项目”包括用户过去操作过的项目，以及足够相似且值得推荐的项目。
初始模型基于按分类体系聚合的订单频率。通过增强模型以模拟随机行为并使用从用户集体参与分类体系中不断学习的分类模型，也解决了各种冷启动问题。
初始模型随后得到增强，以模拟和使用用户偏好，这些偏好可以从用户操作的项目名称中的词汇和其价格中辨别出来。像初始模型一样，用户偏好被映射到分类体系中。这适应了在特定分类体系中的细致学习。例如
-
用户愿意为这家餐厅的主菜支付$8，但不愿为饮料付费。
-
用户喜欢这家餐厅的三文鱼开胃菜，但不喜欢三文鱼主菜。
在这篇文章中，建模和评分是通过 Python 伪代码构建的。在各个阶段，当需要时，伪代码会进行重构。在这篇文章中，Python 伪代码作为沟通和讨论想法的语言。它也可能对那些希望练习实现的读者有用。伪代码的模块化程度足够高，可以逐步实现。
参考文献
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
