开发科学软件
第一部分：测试驱动开发的原则
·
跟随 发表在 Towards Data Science · 10 分钟阅读 · 2023 年 7 月 1 日
–
照片由 Noah Windler提供，拍摄于Unsplash
我们生活在一个计算世界快速扩展可能性的时代。人工智能在解决旧问题和新问题方面不断取得进展，常常以完全意想不到的方式进行。庞大的数据集现在几乎在任何领域都变得普遍，而不仅仅是科学家们在昂贵设施中才能获得的东西。
然而，在过去几十年中，处理数据的软件开发面临的许多挑战依然存在——或者在处理这些新的、大量的数据时问题更加严重。
科学计算领域，传统上专注于开发快速而准确的方法来解决科学问题，近年来已超越其原始的狭窄范围变得相关。在本文中，我将揭示在开发高质量科学软件时出现的一些挑战，以及一些克服这些挑战的策略。我们的最终目标是制定一个逐步指南，以确保准确和高效的开发过程。在后续文章中，我将跟随这个逐步指南解决一个 Python 示例问题。阅读后查看！
TDD 和科学计算：不是天作之合？
测试驱动开发 (TDD)重新定义了软件工程，使开发人员能够编写更耐用、无漏洞的代码。如果你曾经使用过 TDD，你可能对它在编写高质量软件方面的力量非常熟悉。如果没有，希望通过阅读本文你会理解它的重要性。无论你对 TDD 的经验如何，任何熟悉科学计算的人都知道，自动化测试软件的实施往往很棘手。
我推荐大家至少读一次的 TDD 开发周期 提供了一些明智的指示，说明如何以一种每一段代码都通过测试验证的方式开发软件。定期测试可以确保错误常常在被引入之前就被发现。
但 TDD 的一些原则可能与科学软件开发过程完全相悖。例如，在 TDD 中，测试是在代码之前编写的；代码是为了适应测试而编写的。
但设想你正在实现一种全新的数据处理方法。你如何在甚至没有代码的情况下编写测试？TDD 依赖于预期行为：如果在实施新方法之前没有办法量化行为，那么首先编写测试在逻辑上是不可能的！我会认为这种情况很少见，但即使发生了，TDD 仍然可以帮助我们。怎么做呢？
验证与确认
Rilee 和 Clune 观察到（强调部分为我所加）：
有效的数值软件测试需要一套全面的预言器 […] 以及对不可避免的数值误差的稳健估计 […] 初看这些问题常常显得极具挑战性甚至难以克服。然而，我们认为这种普遍看法是不正确的，并由 (1) 模型验证与软件确认的混淆 以及 (2) 科学界倾向于开发相对粗糙的大型程序，聚合了许多算法步骤 驱动。
Oracle 是已知的输入/输出对，这些对可能涉及或不涉及复杂的计算。Oracle 用于传统的 TDD，但它们通常非常简单。它们在科学软件中扮演着更重要的角色，而不仅仅是作为单元测试的一部分！
当我们谈论使用 oracles 检查某些预期行为时，我们谈论的是软件的 验证。对于软件来说，它不关心验证什么，只要输入 X 导致输出 Y。另一方面，确认 是确保代码的输出 Y 准确匹配科学家期望的过程。这个过程 必须必然利用 科学家的领域知识，以实验、模拟、观察、文献调查、数学模型等形式呈现。
这个重要的区别并非科学计算领域所独有。任何 TDD 实践者无论是隐性还是显性地开发了包含验证和确认两个方面的测试。
假设你正在编写代码以将一组人分配到一组标记的椅子上。一个验证测试可以检查 N 个人和 M 把椅子的列表是否输出了 N 个 2 元组。或者检查如果任何列表为空，输出也必须是空列表。与此同时，一个确认测试可以检查如果输入列表包含重复项，函数是否抛出错误。或者检查在任何输出中，没有两个人被分配到同一把椅子上。这些测试需要对我们的问题有领域知识。
虽然 TDD 涉及验证和确认两个方面，但重要的是不要混淆这两者，并在软件开发的适当阶段使用它们。如果你正在编写科学软件——即任何复杂的数值代码，尤其是性能关键的代码——请继续阅读，以了解如何恰当地利用 TDD 达成这些目的。
科学软件中的测试警告
标准软件和科学软件之间的一个重要区别是，在标准软件中，相等性通常是不具争议的。当测试是否有两个人被分配到同一把椅子时，检查标签（可以建模为整数）是否相同是直接的。在科学软件中，浮点数的普遍使用使问题复杂化。相等性不能一般通过 == 检查，通常需要选择数值精度。实际上，精度的定义可能会根据应用而有所不同（例如，见 相对与绝对容差）。这里是一些推荐的数值准确性测试实践：
-
从测试容差开始，精度要尽可能接近所使用的最不精确的浮点类型。你的测试可能会失败。如果失败，请逐位松弛精度，直到测试通过。如果你无法获得良好的精度（例如，你需要
10^-2的容差来使使用 float64 操作的测试通过），可能存在 bug。 -
数值误差通常随着操作次数的增加而增长。在可能的情况下，通过领域特定知识验证精度（例如，泰勒方法具有可以在测试中利用的显式余项，但这种情况很少见）。
-
尽可能偏好绝对容差，避免在比较接近零的值时使用相对容差（“准确性”）。
-
在不同机器上运行测试数千次时，精度单元测试失败并不罕见。如果这种情况持续发生，则要么精度要求过于严格，要么引入了错误。在我的经验中，后一种情况更为常见。
浮点数 😛（照片由 Johannes W 提供，来自 Unsplash）
测试新方法
在开发科学软件时，不能仅仅依赖数值准确性。通常，新方法可以提高准确性或完全改变解决方案，从而从科学家的角度提供“更好的”解决方案。在前一种情况下，科学家可能会使用具有降低容差的先前神谕来确保正确性。在后一种情况下，科学家可能需要创建一个全新的神谕。创建一个精心挑选的神谕示例套件至关重要，这些示例可能会或可能不会被检查数值精度，但科学家可以检查这些示例。
-
精心挑选一组具有代表性的示例，以便可以自动或手动检查。
-
示例应该具有代表性。这可能涉及运行计算密集型任务。因此，重要的是要与单元测试套件解耦。
-
尽可能定期运行这些示例。
随机测试
科学软件可能需要处理非确定性行为。关于如何处理这一点有许多不同的哲学。我个人的方法是通过种子值尽可能控制随机性。这已经成为机器学习实验中的标准，我相信这也是进行通用科学计算的“正确方式”。
我还相信，猴子测试（即模糊测试）——在每次运行时测试随机值的做法——在开发科学软件中扮演着极其宝贵的角色。适当使用猴子测试可以发现隐蔽的错误并增强你的单元测试库。如果使用不当，它可能会创建一个完全不可预测的测试套件。好的猴子测试具有以下特点：
-
测试必须可重现。记录所有重新运行测试所需的种子。
-
随机输入必须覆盖所有可能的输入，并且仅覆盖这些可能的输入。
-
如果可以预测边界情况，则单独处理这些情况。
-
测试应该能够捕捉错误和其他不良行为，除了测试准确性外。如果测试不能标记不良行为，则毫无意义。
-
应该将不良行为作为单独的测试进行研究和隔离，这些测试检查会生成这些错误的整个情况类别（例如，如果输入 -1 失败，并且经过调查发现所有负数都失败，则创建一个针对所有负数的测试）。
分析
除了验证和确认外，开发高性能科学软件的开发人员必须注意性能回归。因此，分析是开发过程的一个重要部分，确保你能从代码中获得最佳性能。
但分析可能很棘手。以下是我在分析科学软件时使用的一些指导原则。
-
分析单元。类似于测试单元，你应该分析性能关键的代码单元。NVIDIA 的 CUDA 最佳实践模型是 评估、并行化、优化、部署 (APOD)。分析单元让你能够评估是否要将代码移植到 GPU 上。
-
首先分析重要的部分。尽量小心，但不要分析不会反复运行的代码片段，或者优化这些片段不会带来显著的收益。
-
多样化分析。分析 CPU 时间、内存和应用程序的任何其他有用指标。
-
确保用于分析的环境是可重复的。包括库版本、CPU 负载等。
-
尝试在单元测试中进行分析。你不需要让回归测试失败，但至少应该标记它们。
将一切整合起来：逐步模型
在这一部分，我们将简要描述我为科学软件应用的开发方法学的主要阶段。这些步骤得到了在学术界、工业界和开源项目中编写科学软件的经验的启发，遵循了上述最佳实践。虽然我不能说我总是应用了这些方法，但我可以诚实地说，我总是后悔没有这样做！
实施周期
-
收集需求。你将如何使用你的方法？考虑它必须提供什么功能、灵活性、输入和输出、是否独立或作为某个更大代码库的一部分。考虑它现在必须做什么以及你将来可能希望它做什么。在这个阶段很容易过早优化，所以记住：“保持简单，傻瓜” 和 “你不会需要它”。
-
勾画设计。创建一个模板，无论是代码还是图示，以建立满足上述要求的设计。
-
实现初始测试。你在第 3 步，迫不及待想开始编码。深呼吸！你会开始编码，但不是你的方法/功能。在这个步骤你编写超简单的测试。真的很小。从简单的验证测试开始，然后转到基本的验证测试。对于验证测试，我建议在开始时尽可能利用分析预设。如果无法做到，请跳过它们。
-
实现你的 alpha 版本。你已经有了测试（验证），可以开始实际实现代码来满足这些测试，而不必担心（很）错误。这个初始实现不必是最快的，但需要是正确的（验证）！我的建议是：从一个简单的实现开始，利用标准库。依赖标准库可以大大降低不正确实现的风险，因为它们利用了自己的测试套件。
-
建立预设库。我不能过分强调这一点！在这一点上，你需要建立可靠的预设，以便你在未来的实现和/或方法变更中始终可以依赖它们。这部分通常在传统的测试驱动开发（TDD）中缺失，但在科学软件中至关重要。它确保你的结果不仅在数值上是正确的，而且可以防止新的或可能不同的实现科学上不准确。在构建验证预设时来回折腾实现和探索脚本是正常的，但避免同时编写测试。
-
重新审视测试。利用你辛勤保存的预设，编写更多的验证单元测试。同样，避免在实现和测试之间来回折腾。
-
实现性能分析。在你的单元测试内部和外部设置性能分析。一旦你完成了第一次迭代，你会回来处理这个问题。
优化周期
-
优化。你现在要使这个函数在你的应用中尽可能快速。利用你的测试和性能分析工具，你可以发挥你的科学计算知识来加速它。
-
重新实现。在这里你可以考虑新的实现方式，例如使用硬件加速库如 GPU、分布式计算等。我建议使用 NVIDIA 的 APOD（评估、并行化、优化、部署）作为一种良好的优化方法论。你可以返回实现周期，但现在你总是拥有一堆预设和测试。如果你预期功能会改变，请参见下面。
新方法周期
-
实现新方法。按照实现周期进行操作，直到第 6 步，包括第 6 步，就像你没有任何预设一样。
-
验证与之前策划的神谕的对比。在神谕构建步骤之后，你可以利用之前实施中的神谕示例，以确保新的神谕在某种程度上“更好”。这一步骤在开发对各种数据都具有鲁棒性的算法和方法中至关重要。在工业界，这一过程被频繁使用，以确保新算法在各种相关情况下表现良好。
下一步
许多这些原则可能只有在查看具体示例时才真正有意义。科学计算涉及多种不同类型的软件，服务于许多目的，因此一种方法很少适用于所有情况。
我鼓励你查看本系列的下一部分，了解如何在实践中实施这些步骤。
开发科学软件
第二部分：Python 的实际应用
·
关注 发表在 Towards Data Science · 14 min read · 2023 年 7 月 1 日
–
在本文中，我们将遵循 本系列第一部分 中提出的 TDD 原则，开发一种称为 Sobel 滤波器的边缘检测滤波器。
在第一篇文章中，我们讨论了为科学软件开发可靠测试方法的重要性——以及它的复杂性。我们还看到如何通过遵循受 TDD 启发但适用于科学计算的开发周期来克服这些问题。下面是这些指令的简化版本。
实现周期
-
收集需求
-
草拟设计
-
实现初步测试
-
实现你的 Alpha 版本
-
构建预言库
-
重新审视测试
-
实现性能分析
优化周期
-
优化
-
重新实现
新方法周期
-
实现新方法
-
与之前策划的预言进行验证
开始：Sobel 滤波器
在本文中，我们将按照上述指令开发一个应用 Sobel 滤波器 的函数。Sobel 滤波器是常用的计算机视觉工具，用于检测或增强图像中的边缘。继续阅读以查看一些示例！
图 1. Sobel–Feldman 算子的核心。来源：自己的工作。
从第一步开始，我们收集一些需求。我们将遵循 这篇文章 中描述的 Sobel 滤波器的标准公式。简单来说，Sobel 算子包括用以下两个 3 × 3 核心对图像 A 进行卷积，平方每个像素，求和并取逐点平方根。如果 Ax 和 Ay 是卷积后的图像，那么 Sobel 滤波器 S 的结果是 √(Ax² + Ay²)。
需求
我们希望这个函数接受任何 2D 数组并生成另一个 2D 数组。我们可能希望它在 ndarray 的任何两个轴上操作。我们甚至可能希望它在两个以上（或以下）轴上工作。我们可能有处理数组边缘的规格。
现在让我们记住保持简单，先从 2D 实现开始。但在此之前，让我们草拟设计。
草拟设计
我们从一个简单的设计开始，利用 Python 的注释。我强烈推荐尽可能多地注释，并使用 mypy 作为 linter。
from typing import Tuple
from numpy.core.multiarray import normalize_axis_index
from numpy.typing import NDArray
def sobel(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:
"""Compute the Sobel filter of an image
Parameters
----------
arr : NDArray
Input image
axes : Tuple[int, int], optional
Axes over which to compute the filter, by default (-2, -1)
Returns
-------
NDArray
Output
"""
# Only accepts 2D arrays
if arr.ndim != 2:
raise NotImplementedError
# Ensure that the axis[0] is the first axis, and axis[1] is the second
# axis. The obscure `normalize_axis_index` converts negative indices to
# indices between 0 and arr.ndim - 1.
if any(
normalize_axis_index(ax, arr.ndim) != i
for i, ax in zip(range(2), axes)
):
raise NotImplementedError
pass
实现测试
这个函数功能不多。但它有文档、注释，并且其当前限制也已包含在内。现在我们有了设计，我们立即转向测试。
首先，我们注意到空图像（全零）没有边缘。因此，它们也必须输出零。实际上，任何常量图像也应返回零。让我们把这点融入到第一个测试中。我们还将探讨如何使用猴子测试来测试多个数组。
# test_zero_constant.py
import numpy as np
import pytest
# Test multiple dtypes at once
@pytest.mark.parametrize(
"dtype",
["float16", "float32", "float64", "float128"],
)
def test_zero(dtype):
# Set random seed
seed = int(np.random.rand() * (2**32 - 1))
np.random.seed(seed)
# Create a 2D array of random shape and fill with zeros
nx, ny = np.random.randint(3, 100, size=(2,))
arr = np.zeros((nx, ny), dtype=dtype)
# Apply sobel function
arr_sob = sobel(arr)
# `assert_array_equal` should fail most of the times.
# It will only work when `arr_sob` is identically zero,
# which is usually not the case.
# DO NOT USE!
# np.testing.assert_array_equal(
# arr_sob, 0.0, err_msg=f"{seed=} {nx=}, {ny=}"
# )
# `assert_almost_equal` can fail when used with high decimals.
# It also relies on float64 checking, which might fail for
# float128 types.
# DO NOT USE!
# np.testing.assert_almost_equal(
# arr_sob,
# np.zeros_like(arr),
# err_msg=f"{seed=} {nx=}, {ny=}",
# decimal=4,
# )
# `assert_allclose` with custom tolerance is my preferred method
# The 10 is arbitrary and depends on the problem. If a method
# which you know to be correct does not pass, increase to 100, etc.
# If the tolerance needed to make the tests pass is too high, make
# sure the method is actually correct.
tol = 10 * np.finfo(arr.dtype).eps
err_msg = f"{seed=} {nx=}, {ny=} {tol=}" # Log seeds and other info
np.testing.assert_allclose(
arr_sob,
np.zeros_like(arr),
err_msg=err_msg,
atol=tol, # rtol is useless for desired=zeros
)
@pytest.mark.parametrize(
"dtype", ["float16", "float32", "float64", "float128"]
)
def test_constant(dtype):
seed = int(np.random.rand() * (2**32 - 1))
np.random.seed(seed)
nx, ny = np.random.randint(3, 100, size=(2,))
constant = np.random.randn(1).item()
arr = np.full((nx, ny), fill_value=constant, dtype=dtype)
arr_sob = sobel(arr)
tol = 10 * np.finfo(arr.dtype).eps
err_msg = f"{seed=} {nx=}, {ny=} {tol=}"
np.testing.assert_allclose(
arr_sob,
np.zeros_like(arr),
err_msg=err_msg,
atol=tol, # rtol is useless for desired=zeros
)
此代码片段可以从命令行运行
$ pytest -qq -s -x -vv --durations=0 test_zero_constant.py
Alpha 版本
当然，我们的测试目前会失败，但现在足够了。让我们实现第一个版本。
from typing import Tuple
import numpy as np
from numpy.core.multiarray import normalize_axis_index
from numpy.typing import NDArray
def sobel(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:
if arr.ndim != 2:
raise NotImplementedError
if any(
normalize_axis_index(ax, arr.ndim) != i
for i, ax in zip(range(2), axes)
):
raise NotImplementedError
# Define our filter kernels. Notice they inherit the input type, so
# that a float32 input never has to be cast to float64 for computation.
# But can you see where using another dtype for Gx and Gy might make
# sense for some input dtypes?
Gx = np.array(
[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]],
dtype=arr.dtype,
)
Gy = np.array(
[[-1, -2, -1], [0, 0, 0], [1, 2, 1]],
dtype=arr.dtype,
)
# Create the output array and fill with zeroes
s = np.zeros_like(arr)
for ix in range(1, arr.shape[0] - 1):
for iy in range(1, arr.shape[1] - 1):
# Pointwise multiplication followed by sum, aka convolution
s1 = (Gx * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
s2 = (Gy * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
s[ix, iy] = np.hypot(s1, s2) # np.sqrt(s1**2 + s2**2)
return s
使用这个新函数，所有测试应该通过，我们应该得到这样的输出：
$ pytest -qq -s -x -vv --durations=0 test_zero_constant.py
........
======================================== slowest durations =========================================
0.09s call t_049988eae7f94139a7067f142bf2852f.py::test_constant[float16]
0.08s call t_049988eae7f94139a7067f142bf2852f.py::test_zero[float64]
0.06s call t_049988eae7f94139a7067f142bf2852f.py::test_constant[float128]
0.04s call t_049988eae7f94139a7067f142bf2852f.py::test_zero[float128]
0.04s call t_049988eae7f94139a7067f142bf2852f.py::test_constant[float64]
0.02s call t_049988eae7f94139a7067f142bf2852f.py::test_constant[float32]
0.01s call t_049988eae7f94139a7067f142bf2852f.py::test_zero[float16]
(17 durations < 0.005s hidden. Use -vv to show these durations.)
8 passed in 0.35s
到目前为止，我们有：
-
收集了我们问题的需求。
-
绘制了初步设计。
-
实现了一些测试。
-
实现了通过这些测试的 alpha 版本。
这些测试为未来版本提供了验证，以及一个非常基础的oracle 库。但是，尽管单元测试在捕捉与预期结果的细微偏差方面非常出色，它们在区分错误结果与数量上不同——但仍然正确——的结果方面表现并不好。假设明天我们要实现另一个 Sobel 型操作符，它有一个更长的内核。我们不期望结果与当前版本完全匹配，但我们期望两个函数的输出在质量上非常相似。
此外，尝试对我们的函数进行多种不同输入是一个很好的做法，以了解它们对不同输入的表现。这确保了我们科学地验证结果。
Scikit-image拥有一个出色的图像库，我们可以用它来创建 oracles。
# !pip installscikit-image pooch
from typing import Dict, Callable
import numpy as np
import skimage.data
bwimages: Dict[str, np.ndarray] = {}
for attrname in skimage.data.__all__:
attr = getattr(skimage.data, attrname)
# Data are obtained via function calls
if isinstance(attr, Callable):
try:
# Download the data
data = attr()
# Ensure it is a 2D array
if isinstance(data, np.ndarray) and data.ndim == 2:
# Convert from various int types to float32 to better
# assess precision
bwimages[attrname] = data.astype(np.float32)
except:
continue
# Apply sobel to images
bwimages_sobel = {k: sobel(v) for k, v in bwimages.items()}
一旦我们有了数据，就可以绘制它。
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
def create_colorbar(im, ax):
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cb = ax.get_figure().colorbar(im, cax=cax, orientation="vertical")
return cax, cb
for name, data in bwimages.items():
fig, axs = plt.subplots(
1, 2, figsize=(10, 4), sharex=True, sharey=True
)
im = axs[0].imshow(data, aspect="equal", cmap="gray")
create_colorbar(im, axs[0])
axs[0].set(title=name)
im = axs[1].imshow(bwimages_sobel[name], aspect="equal", cmap="gray")
create_colorbar(im, axs[1])
axs[1].set(title=f"{name} sobel")
图 2. “Binary blobs”数据集在 Sobel 滤波之前（左）和之后（右）。来源：自制。
图 3. “Text”数据集在 Sobel 滤波之前（左）和之后（右）。来源：自制。
这些看起来非常好！我建议将这些数据存储为数组和图形，以便我能快速对照新版本。我强烈推荐HD5F用于数组存储。你还可以设置一个Sphinx Gallery，在更新文档时直接生成图形，这样你就有了一个可重复的图形构建，可以用来与以前的版本进行对比。
在结果经过验证后，你可以将它们存储在磁盘上，并将它们作为单元测试的一部分进行使用。类似这样：
oracle_library = [(k, v, bwimages_sobel[k]) for k, v in bwimages.items()]
# save_to_disk(oracle_library, ...)
# test_oracle.py
import numpy as np
import pytest
from numpy.typing import NDArray
# oracle_library = read_from_disk(...)
@pytest.mark.parametrize("name,input,output", oracle_library)
def test_oracles(name: str, input: NDArray, output: NDArray):
output_new = sobel(input)
tol = 10 * np.finfo(input.dtype).eps
mean_avg_error = np.abs(output_new - output).mean()
np.testing.assert_allclose(
output_new,
output,
err_msg=f"{name=} {tol=} {mean_avg_error=}",
atol=tol,
rtol=tol,
)
性能分析
计算这些数据集的 Sobel 滤波器花费了一段时间！因此下一步是对代码进行性能分析。我建议为每个测试创建一个benchmark_xyz.py文件，并定期重新运行它们，以探测性能回归。这甚至可以成为你的单元测试的一部分，但在这个示例中我们不会深入探讨。另一个想法是使用上述 oracles 进行基准测试。
有许多方法来计时代码执行。要获得系统范围的“真实”经过时间，使用内置time模块中的perf_counter_ns来以纳秒为单位测量时间。在 Jupyter notebook 中，内置的[%%timeit](https://ipython.readthedocs.io/en/stable/interactive/magics.html#cell-magics) cell magic会对某个单元执行时间进行计时。下面的装饰器灵感来源于这个 cell magic，并可以用来计时任何函数。
import time
from functools import wraps
from typing import Callable, Optional
def sizeof_fmt(num, suffix="s"):
for unit in ["n", "μ", "m"]:
if abs(num) < 1000:
return f"{num:3.1f} {unit}{suffix}"
num /= 1000
return f"{num:.1f}{suffix}"
def timeit(
func_or_number: Optional[Callable] = None,
number: int = 10,
) -> Callable:
"""Apply to a function to time its executions.
Parameters
----------
func_or_number : Optional[Callable], optional
Function to be decorated or `number` argument (see below), by
default None
number : int, optional
Number of times the function will run to obtain statistics, by
default 10
Returns
-------
Callable
When fed a function, returns the decorated function. Otherwise
returns a decorator.
Examples
--------
.. code-block:: python
@timeit
def my_fun():
pass
@timeit(100)
def my_fun():
pass
@timeit(number=3)
def my_fun():
pass
"""
if isinstance(func_or_number, Callable):
func = func_or_number
number = number
elif isinstance(func_or_number, int):
func = None
number = func_or_number
else:
func = None
number = number
def decorator(f):
@wraps(f)
def wrapper(*args, **kwargs):
runs_ns = np.empty((number,))
# Run first and measure store the result
start_time = time.perf_counter_ns()
result = f(*args, **kwargs)
runs_ns[0] = time.perf_counter_ns() - start_time
for i in range(1, number):
start_time = time.perf_counter_ns()
f(*args, **kwargs) # Without storage, faster
runs_ns[i] = time.perf_counter_ns() - start_time
time_msg = f"{sizeof_fmt(runs_ns.mean())} ± "
time_msg += f"{sizeof_fmt(runs_ns.std())}"
print(
f"{time_msg} per run (mean ± std. dev. of {number} runs)"
)
return result
return wrapper
if func is not None:
return decorator(func)
return decorator
对我们的函数进行测试：
arr_test = np.random.randn(500, 500)
sobel_timed = timeit(sobel)
sobel_timed(arr_test);
# 3.9s ± 848.9 ms per run (mean ± std. dev. of 10 runs)
不太快 😦
在调查慢速或性能回归时，还可以跟踪每一行的运行时间。[line_profiler](https://github.com/pyutils/line_profiler) 库 是一个很好的资源。它可以通过 Jupyter 单元魔法 或使用 API 来使用。以下是一个 API 示例：
from line_profiler import LineProfiler
lp = LineProfiler()
lp_wrapper = lp(sobel)
lp_wrapper(arr_test)
lp.print_stats(output_unit=1) # 1 for seconds, 1e-3 for milliseconds, etc.
这是一个示例输出：
Timer unit: 1 s
Total time: 4.27197 s
File: /tmp/ipykernel_521529/1313985340.py
Function: sobel at line 8
Line # Hits Time Per Hit % Time Line Contents
==============================================================
8 def sobel(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:
9 # Only accepts 2D arrays
10 1 0.0 0.0 0.0 if arr.ndim != 2:
11 raise NotImplementedError
12
13 # Ensure that the axis[0] is the first axis, and axis[1] is the second
14 # axis. The obscure `normalize_axis_index` converts negative indices to
15 # indices between 0 and arr.ndim - 1.
16 1 0.0 0.0 0.0 if any(
17 normalize_axis_index(ax, arr.ndim) != i
18 1 0.0 0.0 0.0 for i, ax in zip(range(2), axes)
19 ):
20 raise NotImplementedError
21
22 1 0.0 0.0 0.0 Gx = np.array(
23 1 0.0 0.0 0.0 [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]],
24 1 0.0 0.0 0.0 dtype=arr.dtype,
25 )
26 1 0.0 0.0 0.0 Gy = np.array(
27 1 0.0 0.0 0.0 [[-1, -2, -1], [0, 0, 0], [1, 2, 1]],
28 1 0.0 0.0 0.0 dtype=arr.dtype,
29 )
30 1 0.0 0.0 0.0 s = np.zeros_like(arr)
31 498 0.0 0.0 0.0 for ix in range(1, arr.shape[0] - 1):
32 248004 0.1 0.0 2.2 for iy in range(1, arr.shape[1] - 1):
33 248004 1.8 0.0 41.5 s1 = (Gx * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
34 248004 1.7 0.0 39.5 s2 = (Gy * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
35 248004 0.7 0.0 16.8 s[ix, iy] = np.hypot(s1, s2)
36 1 0.0 0.0 0.0 return s
大量时间花费在最内层的循环中。NumPy 更喜欢矢量化数学，因为这样可以依赖编译的代码。由于我们使用显式的 for 循环，因此这个函数非常慢是合理的。
此外，明智地处理循环中的内存分配非常重要。NumPy 对于在循环中分配小量内存有一定的智能，因此定义s1或s2的行可能并未分配新的内存。但它们也可能分配了，或者可能在幕后有一些我们未察觉的内存分配。因此，我建议也进行内存分析。我喜欢使用 Memray，但还有其他工具，如 Fil 和 Sciagraph。我还会避免使用 memory_profiler，因为（非常不幸的是！）它不再维护。此外，Memray 更强大。以下是 Memray 使用中的一个示例：
$ # Create sobel.bin which holds the profiling information
$ memray run -fo sobel.bin --trace-python-allocators sobel_run.py
Writing profile results into sobel.bin
Memray WARNING: Correcting symbol for aligned_alloc from 0x7fc5c984d8f0 to 0x7fc5ca4a5ce0
[memray] Successfully generated profile results.
You can now generate reports from the stored allocation records.
Some example commands to generate reports:
python3 -m memray flamegraph sobel.bin
$ # Generate flame graph
$ memray flamegraph -fo sobel_flamegraph.html --temporary-allocations sobel.bin
⠙ Calculating high watermark... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 99% 0:00:0100:01
⠏ Processing allocation records... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 98% 0:00:0100:01
Wrote sobel_flamegraph.html
$ # Show memory tree
$ memray tree --temporary-allocations sobel.bin
⠧ Calculating high watermark... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 100% 0:00:0100:01
⠧ Processing allocation records... ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 100% 0:00:0100:01
Allocation metadata
-------------------
Command line arguments:
'memray run -fo sobel.bin --trace-python-allocators sobel_run.py'
Peak memory size: 11.719MB
Number of allocations: 15332714
Biggest 10 allocations:
-----------------------
📂 123.755MB (100.00 %) <ROOT>
└── [[3 frames hidden in 2 file(s)]]
└── 📂 123.755MB (100.00 %) _run_code /usr/lib/python3.10/runpy.py:86
├── 📂 122.988MB (99.38 %) <module> sobel_run.py:40
│ ├── 📄 51.087MB (41.28 %) sobel sobel_run.py:35
│ ├── [[1 frames hidden in 1 file(s)]]
│ │ └── 📄 18.922MB (15.29 %) _sum
│ │ lib/python3.10/site-packages/numpy/core/_methods.py:49
│ └── [[1 frames hidden in 1 file(s)]]
│ └── 📄 18.921MB (15.29 %) _sum
│ lib/python3.10/site-packages/numpy/core/_methods.py:49
...
图 4. Memray 火焰图（alpha 版本）。来源：自制作品。
Beta 版本和一个 Bug
现在我们有了一个可用的 alpha 版本和一些性能分析函数，我们将利用 SciPy 库 来获得一个更快的版本。
from typing import Tuple
import numpy as np
from numpy.core.multiarray import normalize_axis_index
from numpy.typing import NDArray
from scipy.signal import convolve2d
def sobel_conv2d(
arr: NDArray, axes: Tuple[int, int] = (-2, -1)
) -> NDArray:
if arr.ndim != 2:
raise NotImplementedError
if any(
normalize_axis_index(ax, arr.ndim) != i
for i, ax in zip(range(2), axes)
):
raise NotImplementedError
# Create the kernels as a single, complex array. Allows us to use
# np.abs instead of np.hypot to calculate the magnitude.
G = np.array(
[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]],
dtype=arr.dtype,
)
G = G + 1j * np.array(
[[-1, -2, -1], [0, 0, 0], [1, 2, 1]],
dtype=arr.dtype,
)
s = convolve2d(arr, G, mode="same")
np.absolute(s, out=s) # In-place abs
return s.real
sobel_timed = timeit(sobel_conv2d)
sobel_timed(arr_test)
# 14.3 ms ± 1.71 ms per loop (mean ± std. dev. of 10 runs)
好得多！但它是正确的吗？
图 5. “微动脉瘤” 数据集在使用两个版本的 Sobel 滤波前（左）和后（中间和右）。来源：自制作品。
图像看起来非常相似，但如果注意颜色比例，它们并不相同。运行测试时标记了一个小的平均误差。幸运的是，我们现在可以很好地检测定量和定性差异。
在调查这个 Bug 后，我们将其归因于不同的边界条件。查看 [convolve2d](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.convolve2d.html) 文档 告诉我们，输入数组在应用卷积核之前会用零填充。在 alpha 版本中，我们填充了输出！
哪一个是正确的？可以说 SciPy 实现更有意义。在这种情况下，我们应该采用新版本作为标准，如果需要，修复 alpha 版本以匹配它。这在科学软件开发中很常见：如何更好地做事的新信息会改变标准和测试。
在这种情况下，修复很简单，在处理之前用零填充数组即可。
def sobel_v2(arr: NDArray, axes: Tuple[int, int] = (-2, -1)) -> NDArray:
# ...
arr = np.pad(arr, (1,)) # After padding, it is shaped (nx + 2, ny + 2)
s = np.zeros_like(arr)
for ix in range(1, arr.shape[0] - 1):
for iy in range(1, arr.shape[1] - 1):
s1 = (Gx * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
s2 = (Gy * arr[ix - 1 : ix + 2, iy - 1 : iy + 2]).sum()
s[ix - 1, iy - 1] = np.hypot(s1, s2) # Adjust indices
return s
一旦我们修正了函数，就可以更新依赖于它们的测试和预期结果。
最后的思考
我们展示了如何在这篇文章中探讨的一些软件开发理念中付诸实践。我们还介绍了一些工具，你可以使用它们来开发高质量、高性能的代码。
我建议你在自己的项目中尝试这些想法。特别是，练习代码分析和改进。我们编写的 Sobel 滤波器函数非常低效，我建议尝试改进它。
例如，可以尝试使用如Numba的 JIT 编译器进行 CPU 并行化，将内部循环移植到Cython，或者使用Numba或CuPy实现 CUDA GPU 函数。欢迎查看我关于用 Numba 编写 CUDA 内核的系列文章。
细节决定成败：通过跳出框框思维成为 Power BI 大师
Power BI 中充满了“无名英雄”！其中之一，分析面板，结合更改可视化类型，帮助我显著提升了 Power BI 报告的性能
·发布于 Towards Data Science ·阅读时间 5 分钟·2023 年 7 月 7 日
–
照片由 Alice Dietrich 提供，来源于 Unsplash
几周前，我正在为我的一个客户进行Power BI 报告的性能调优。报告页面的渲染非常缓慢（超过 15 秒）。为了给你一点背景：该报告使用了与托管在 SSAS Tabular 2016 中的表格模型的实时连接。
如果我告诉你，我在没有更改任何 DAX 代码的情况下，将报告页面的性能提升了两倍，你会如何反应？
继续阅读，你会发现为什么细节往往隐藏着问题，以及跳出框框思维如何帮助你成为真正的 Power BI 大师:)
图片来源：作者
让我快速解释上面的插图。这里有一个折线图和簇状柱形图，其中四条线代表用户从左侧的切片器中选择（报告阈值和三个层次），而柱形图则是总销售额。数据按年份和产品进行细分。每条线都使用 DAX 计算（顺便提一下，SSAS 2016 中没有SELECTEDVALUE 函数）。
我已启用性能分析器，获取了 DAX 查询并在 DAX Studio 中执行了它：
图片来源：作者
如你所见，查询执行时间为 13.5 秒（在运行前清除缓存），而大部分时间花费在公式引擎中（76%）。这很重要，因为我们将把这个结果与改进后的报告页面进行比较。
那么，一个有经验的 Power BI 开发人员会怎么做来优化这个场景？重写 DAX？错误！
让我们看看在不更改 DAX 逻辑的情况下可以做些什么！
如果你不知道，Power BI 提供了一个被低估的功能，或者我们可以称之为“无名英雄”，即分析面板。
作者提供的图片
显然，我们中的大多数人将 Power BI 开发时间花在了另外两个面板——数据和格式上。所以，你会惊讶于有多少 Power BI 开发人员甚至不知道这个第三个面板，或者即使知道，也很少使用它。
不深入细节，这个面板让你可以向你的可视化中添加额外的分析成分——如最小值、最大值、平均值、中位数线、误差条等。根据可视化的类型，并非所有选项都可用！在我的使用案例中，这一点非常重要。
当我打开分析面板时，只能看到误差条：
作者提供的图片
本质上，这里的想法是，由于这四条线的值不会根据可视化中的数字变化（它们的值是基于切片器选择的常量值），因此利用分析面板中的常量线功能。由于在折线图和聚合柱状图中没有常量线选项，我们将复制我们的可视化并将其类型更改为常规的聚合柱状图。
作者提供的图片
如你所见，“数字”在这里，但我们的线条却不见了。让我们切换到分析面板，创建 4 条常量线，每条基于切片器选择生成的 DAX 度量：
作者提供的图片
第一步是添加一条常量线。接下来，展开“折线”属性，并选择“fX”按钮，允许你基于表达式（在我们这里是由 DAX 度量生成的表达式）设置常量线的值。对所有四条线重复此过程。
我再提醒你一次，请注意我完全没有触及 DAX 代码！
一旦我关闭了 Y 轴，这就是我的“副本”可视化的样子：
作者提供的图片
基本上是一样的，对吧？
好的，现在让我们检查一下这个可视化的性能：
作者提供的图片
比原始版本快了超过 5 秒！如果你仔细对比之前的 DAX Studio 截图，你会注意到 存储引擎查询 的数量与之前的情况完全相同（而 SE 时间几乎相同），这意味着存储引擎需要完成的工作量完全相同以检索数据。
关键区别在于公式引擎的时间——与原始可视化中 75% 的总查询时间花费在 FE 上不同，这次减少到不到 60%！
我很好奇为什么会发生这种情况，以及公式引擎生成的两个查询计划之间的主要区别是什么。
较慢查询——DAX 代码部分
较快查询——DAX 代码部分
两个查询计划之间的“唯一”区别——在较慢的版本中，创建了两个虚拟表：一个用于计算可视化中的“列”值（_ScopedCoreI0），另一个用于计算同一可视化中的行值（_ScopedCoreDM0）。最后，这两个表使用 NATURALLEFTOUTERJOIN 函数连接。
在更快的版本中，没有计算行值的第二个表。此外，计算行值的度量值被包装在 IGNORE 函数中，该函数标记 SUMMARIZECOLUMNS 表达式中的度量值以便从非空行的评估中省略。
结论
正如你所见，更改可视化类型，加上“无名英雄”分析面板的使用，确保了在这种情况下显著的性能提升。人们常说“魔鬼在细节中”绝非偶然——因此，往往跳出框框思考会带来一些创新的解决方案。
感谢阅读！
训练 LLMs 的不同方式
为什么提示不属于其中任何一种
·
关注 发布于 Towards Data Science ·10 分钟阅读·2023 年 7 月 21 日
–
在大型语言模型（LLMs）领域中，存在多种训练机制，具有不同的方式、要求和目标。由于它们服务于不同的目的，因此重要的是不要混淆它们，并了解它们适用的不同场景。
在本文中，我想概述一些最重要的训练机制，包括预训练、微调、基于人类反馈的强化学习（RLHF）和适配器。此外，我还将讨论提示的作用，尽管它本身不被视为学习机制，并阐明提示调优的概念，它在提示与实际训练之间架起了一座桥梁。
预训练
“训练”。就像“训练”一样。我很抱歉……照片由 Brian Suman 拍摄，来源于 Unsplash
预训练是最基础的训练方式，相当于你可能知道的其他机器学习领域的训练。在这里，你从一个未经训练的模型（即权重随机初始化的模型）开始，并训练模型在给定一系列前序标记的情况下预测下一个标记。为此，从各种来源收集大量句子，并将其分成小块输入模型。
这里使用的训练模式称为 自监督。从被训练模型的角度来看，我们可以说这是一个监督学习方法，因为模型在做出预测后总是能得到正确的答案。例如，给定序列 I like ice … 模型可能预测 cones 作为下一个词，然后被告知答案是错误的，因为实际的下一个词是 cream。最终，可以计算损失，并调整模型权重，以便下次预测得更好。之所以称之为 自 监督（而不是简单地称为 监督），是因为不需要提前通过昂贵的程序来收集标签，而这些标签已经包含在数据中。给定句子 I like ice cream，我们可以自动将其分割为 I like ice 作为输入，cream 作为标签，这不需要人工干预。尽管这不是模型本身完成的，但仍由机器自动执行，因此在学习过程中，AI 自行监督的概念依然存在。
最终，通过大量文本的训练，模型学习了编码语言结构的一般规则（例如，它学习到 I like 可以被名词或分词跟随）以及文本中包含的知识。例如，它学习到，句子 Joe Biden is … 通常会跟随 the president of the United States，因此对这一知识有了表征。
这种预训练已经由其他人完成，你可以直接使用像 GPT 这样的模型。那么，为什么你还需要训练一个类似的模型呢？如果你处理的数据具有类似语言的特性，但不是普通语言本身，训练一个从头开始的模型可能是必要的。音乐记谱法可能是一个例子，它在某种程度上像语言一样有结构。存在某些规则和模式，关于哪些乐段可以跟随其他乐段，但一个训练在自然语言上的大型语言模型（LLM）不能处理这种数据，所以你需要训练一个新模型。然而，LLM 的架构可能适用，因为音乐记谱法和自然语言之间有许多相似之处。
微调
这些旋钮用于微调弦乐器。照片由 Tony Woodhead 提供，来源于 Unsplash
尽管一个预训练的 LLM 能够执行各种任务，这主要是由于它编码的知识，但它也存在两个主要的缺陷：输出结构和数据中缺失的知识。
正如你所知道的，LLM 总是根据前面的令牌序列预测下一个令牌。对于继续一个给定的故事，这可能没问题，但在其他情况下，这可能不是你想要的。如果你需要不同的输出结构，有两种主要方法可以实现这一目标。你可以通过设计提示，使得模型预测下一个令牌的固有能力完成你的任务（这称为提示工程），或者你可以改变最后一层的输出，使其反映你的任务，就像你在任何其他机器学习模型中一样。比如说分类任务，其中有N个类别。通过提示工程，你可以指示模型在给定输入后总是输出分类标签。通过微调，你可以将最后一层修改为N个输出神经元，并从激活度最高的神经元中得出预测的类别。
LLM 的另一个局限在于它所训练的数据。由于数据来源相当丰富，最著名的 LLM 编码了大量的常识。因此，它们可以告诉你，比如说，美国总统、贝多芬的主要作品、量子物理学的基本原理以及弗洛伊德的主要理论。然而，也有一些领域模型并不了解，如果你需要处理这些领域，微调可能对你有帮助。
微调的思想是利用一个已经预训练的模型，继续用不同的数据进行训练，并且在训练过程中只调整最后几层的权重。这仅需要初始训练所需资源的一小部分，因此可以更快地完成。另一方面，模型在预训练过程中学到的结构仍然编码在前几层中，并且可以加以利用。比如说，你想要教你的模型关于你喜欢的但不太知名的奇幻小说，这些小说并未包含在训练数据中。通过微调，你可以利用模型对自然语言的一般知识，让它理解新的奇幻小说领域。
RLHF 微调
RLHF 微调是关于最大化奖励的。照片由 Alexander Grey 提供，来源于 Unsplash
微调模型的一个特殊情况是基于人类反馈的强化学习（RLHF），这是 GPT 模型和像 Chat-GPT 这样的聊天机器人之间的主要区别之一。通过这种微调，模型被训练以产生人类在与模型对话时认为最有用的输出。
主要思想如下：给定一个任意提示，从模型中生成多个输出。人类根据这些输出的有用性或适当性对其进行排名。给定四个样本 A、B、C 和 D，人类可能会决定 C 是最佳输出，B 略差但与 D 相等，而 A 是最差的输出。这将导致排序 C > B = D > A。接下来，这些数据用于训练一个奖励模型。这是一个全新的模型，它通过给予奖励来学习对 LLM 的输出进行评分，这些奖励反映了人类的偏好。一旦奖励模型训练完成，它就可以代替人类进行评分。现在，模型的输出由奖励模型进行评分，这个奖励作为反馈提供给 LLM，然后 LLM 会根据最大化奖励进行调整；这个想法与 GANs 非常相似。
正如你所见，这种训练需要人工标注的数据，这需要相当大的努力。然而，所需的数据量是有限的，因为奖励模型的思想是从这些数据中进行概括，使其能够在学习其部分后独立对 llm 进行评分。RLHF 常用于使 LLM 的输出更具对话性，或避免模型表现出不良行为，如模型变得刻薄、侵入或侮辱。
适配器
两种适配器可以插入到已经存在的网络中。图片来自 arxiv.org/pdf/2304.01933.pdf。
在上述的微调过程中，我们调整了模型最后几层的一些参数，而之前几层的其他参数保持不变。然而，还有一种替代方案，通过减少所需训练的参数数量来提高效率，这就是所谓的adapters*。
使用适配器意味着在已经训练好的模型上添加额外的层。在微调过程中，只有这些适配器被训练，而模型的其余参数完全不变。然而，这些层比模型自带的层要小得多，这使得调整它们变得更容易。此外，它们可以插入模型的不同位置，而不仅仅是最后的位置。在上面的图片中，你可以看到两个示例：一个是将适配器作为串行层添加，另一个是将其并行添加到已经存在的层中。
提示
提示更多的是告诉模型要做什么，而不是如何去做。照片由 Ian Taylor 提供，发布在 Unsplash
你可能会想知道提示是否算作另一种训练模型的方法。提示是构建实际模型输入之前的指令，特别是如果你使用少量示例提示，你在提示中提供给大语言模型的示例，这与训练非常相似，训练也包含展示给模型的示例。然而，提示与训练模型有所不同。首先，从定义上讲，我们只在更新权重时才谈论训练，而提示过程中不会更新权重。创建提示时，你不会改变任何模型，也不会改变权重，不会产生新的模型，也不会改变模型中编码的知识或表示。提示应该被视为一种指导大语言模型的方法，告诉它你期望从它那里得到什么。以下提示作为一个例子：
"""Classify a given text regarding its sentiment.
Text: I like ice cream.
Sentiment: negative
Text: I really hate the new AirPods.
Sentiment: positive
Text: Donald is the biggest jerk on earth. I hate him so much!
Sentiment: neutral
Text: {user_input}
Sentiment:"""
我指示模型进行情感分类，正如你可能已经注意到的，我给模型的示例全都是错误的！如果一个模型用这样的数据进行训练，它会混淆标签 positive（正面）、negative（负面）和 neutral（中性）。现在，如果我让模型对句子 I like ice cream 进行分类，而这个句子是我给出的示例之一，会发生什么？有趣的是，它将其分类为 positive（正面），这与提示相反，但在语义上是正确的。这是因为提示没有训练模型，也没有改变它所学到的表示。提示只是告诉模型我期望的结构，即我期望情感标签（可以是 positive、negative 或 neutral）跟在冒号后面。
提示调优
提示调优也叫做软提示。软如美利奴羊毛……照片由 Advocator SY 提供，发布在 Unsplash
尽管提示本身并没有训练大语言模型，但有一种机制叫做 提示调优（也叫 软提示），它与提示相关，可以视为一种训练。
在前面的例子中，我们将提示视为一种自然语言文本，这种文本提供给模型以告诉它要做什么，并且在实际输入之前。这就是说，模型输入变成了，例如，<将以下内容标记为正面、负面或中性：> <我喜欢冰淇淋>。当我们自己创建提示时，我们称之为硬提示*。在软提示中，格式将保持不变，但提示本身不是由我们设计的，而是通过数据学习的。具体来说，提示由一个向量空间中的参数组成，这些参数可以在训练过程中进行调优，以获得更小的损失，从而得到更好的回答。也就是说，在训练之后，提示将是那种字符序列，这种序列能为我们提供的数据得到最佳回答。然而，模型参数本身没有经过训练。
提示调优的一个大优点是，你可以为不同的任务训练多个提示，但仍然可以在相同的模型中使用它们。就像在硬提示中，你可能会构造一个用于文本摘要的提示，一个用于情感分析的提示，以及一个用于文本分类的提示，但在同一个模型中使用所有这些提示，你可以为这些目的调优三个提示，并且仍然使用相同的模型。如果你使用微调，恰恰相反，你最终会得到三个只服务于特定任务的模型。
总结
我们刚刚看到了各种不同的训练机制，因此让我们在最后做一个简短的总结。
-
预训练 LLM 意味着以自我监督的方式教它预测下一个标记。
-
微调是调整预训练 LLM 最后几层的权重，可以用来将模型适应特定的上下文。
-
RLHF 的目标是调整模型的行为以匹配人类的期望，并且需要额外的标注工作。
-
适配器通过向预训练的 LLM 添加小层，提供了一种更高效的微调方式。
-
提示不被认为是训练，因为它不会改变模型的内部表示。
-
提示调优是一种调整生成提示的权重的技术，但不会影响模型的权重本身。
当然，还有许多其他的训练机制，每天都有新的机制被发明。LLM 能做的远不止预测文本，教它们做到这些需要多种技能和技术，其中一些我刚刚向你介绍过。
进一步阅读
Instruct-GPT 是 RLHF 最著名的例子之一：
-
Ouyang, L., Wu, J., Jiang, X., Almeida, D., Wainwright, C., Mishkin, P., … & Lowe, R. (2022). 训练语言模型以根据人类反馈遵循指令。神经信息处理系统进展, 35, 27730–27744.
常见适配器形式的概述可以在 LLM-Adapters 项目中找到：
-
Hu, Z., Lan, Y., Wang, L., Xu, W., Lim, E. P., Lee, R. K. W., … & Poria, S. (2023). LLM-Adapters: 一种用于大规模语言模型参数高效微调的适配器家族。 arXiv 预印本 arXiv:2304.01933。
提示调优的探索可以在这里找到：
- Lester, B., Al-Rfou, R., & Constant, N. (2021). 规模的力量在参数高效提示调优中的作用。 arXiv 预印本 arXiv:2104.08691。
一些关于提示调优的不错解释可以在这里找到：
喜欢这篇文章吗？ 关注我 以便获取我未来的帖子通知。
作为 Pytorch 神经网络层的微分方程
原文：
towardsdatascience.com/differential-equations-as-a-pytorch-neural-network-layer-7614ba6d587f
如何在 pytorch 中使用微分方程层
·发表于 Towards Data Science ·阅读时间 17 分钟·2023 年 4 月 26 日
–
微分方程是现代科学的大多数数学基础。它们通过描述变化率（微分）的方程来描述系统的状态。许多系统都可以通过这种形式的方程很好地描述。例如，描述运动、电磁学和量子力学的物理定律都采用这种形式。更广泛地说，微分方程通过质量作用定律描述化学反应速率，通过 SIR 模型描述神经元放电和疾病传播。
深度学习革命带来了新一套工具，用于在巨大的数据集上进行大规模优化。在这篇文章中，我们将探讨如何使用这些工具来拟合 pytorch 中自定义微分方程层的参数。
我们要解决的问题是什么？
假设我们有一些时间序列数据 y(t)，我们希望用微分方程对其建模。数据表现为在时间 tᵢ 上的一组观察值 yᵢ。基于对基础系统的某些领域知识，我们可以写出一个微分方程来近似该系统。
在最一般的形式下，这表现为：
一般常微分方程系统
其中 y 是系统的状态，t 是时间，而 θ 是模型的参数。在这篇文章中，我们将假设参数 θ 是未知的，我们希望从数据中学习这些参数。
本文的所有代码可在 github 或 colab notebook 上找到，所以如果你想跟随学习，无需尝试复制粘贴。
让我们导入本帖所需的库。我们将使用的唯一非标准机器学习库是 torchdiffeq 库来解决微分方程。该库在 pytorch 中实现了数值微分方程求解器。
import torch
import torch.nn as nn
from torchdiffeq import odeint as odeint
import pylab as plt
from torch.utils.data import Dataset, DataLoader
from typing import Callable, List, Tuple, Union, Optional
from pathlib import Path
模型
我们建模过程的第一步是定义模型。对于微分方程，这意味着我们必须选择函数 f(y,t;θ) 的形式和表示参数 θ 的方式。我们还需要以与 pytorch 兼容的方式完成这项工作。
这意味着我们需要将我们的函数编码为 torch.nn.Module 类。正如你所看到的，这非常简单，只需要定义两个方法。让我们从三个示例模型中的第一个开始。
van Der Pol 振荡器 (VDP)
我们可以使用 torch.nn.Module 类来定义一个微分方程系统，其中参数通过 torch.nn.Parameter 声明创建。这让 pytorch 知道我们希望为这些参数累积梯度。我们还可以通过不使用此声明来包含固定参数（我们不想调整的参数）。
我们将使用的第一个示例是经典的 VDP 振荡器，它是一个具有单个参数 μ 的非线性振荡器。该系统的微分方程是：
VDP 振荡器方程
其中 x 和 y 是状态变量。VDP 模型用于模拟从电子电路到心律失常和昼夜节律的一切。我们可以在 pytorch 中定义此系统如下：
class VDP(nn.Module):
"""
Define the Van der Pol oscillator as a PyTorch module.
"""
def __init__(self,
mu: float, # Stiffness parameter of the VDP oscillator
):
super().__init__()
self.mu = torch.nn.Parameter(torch.tensor(mu)) # make mu a learnable parameter
def forward(self,
t: float, # time index
state: torch.TensorType, # state of the system first dimension is the batch size
) -> torch.Tensor: # return the derivative of the state
"""
Define the right hand side of the VDP oscillator.
"""
x = state[..., 0] # first dimension is the batch size
y = state[..., 1]
dX = self.mu*(x-1/3*x**3 - y)
dY = 1/self.mu*x
# trick to make sure our return value has the same shape as the input
dfunc = torch.zeros_like(state)
dfunc[..., 0] = dX
dfunc[..., 1] = dY
return dfunc
def __repr__(self):
"""Print the parameters of the model."""
return f" mu: {self.mu.item()}"
你只需要定义 init 方法 (init) 和 forward 方法。我添加了一个字符串方法 *repr* 来美观地打印参数。关键点在于我们如何在 forward 方法中将微分方程转换为 torch 代码。此方法需要定义微分方程的右侧。
让我们看看如何使用来自 torchdiffeq 的 odeint 方法来集成这个模型：
vdp_model = VDP(mu=0.5)
# Create a time vector, this is the time axis of the ODE
ts = torch.linspace(0,30.0,1000)
# Create a batch of initial conditions
batch_size = 30
# Creates some random initial conditions
initial_conditions = torch.tensor([0.01, 0.01]) + 0.2*torch.randn((batch_size,2))
# Solve the ODE, odeint comes from torchdiffeq
sol = odeint(vdp_model, initial_conditions, ts, method='dopri5').detach().numpy()
plt.plot(ts, sol[:,:,0], lw=0.5);
plt.title("Time series of the VDP oscillator");
plt.xlabel("time");
plt.ylabel("x");p
这里是解的相平面图（动态状态的参数图的相平面图）。
# Check the solution
plt.plot(sol[:,:,0], sol[:,:,1], lw=0.5);
plt.title("Phase plot of the VDP oscillator");
plt.xlabel("x");
plt.ylabel("y");
颜色表示我们批次中的 30 条独立轨迹。解返回为具有 (time_points, batch number, dynamical_dimension) 维度的 torch 张量。
sol.shape
(1000, 30, 2)
Lotka-Volterra 捕食者-猎物方程
作为另一个示例，我们创建一个 Lotka-Volterra 捕食者-猎物方程的模块。在 Lotka-Volterra (LV) 捕食者-猎物模型中，有两个主要变量：猎物的种群 (x) 和捕食者的种群 (y)。该模型由以下方程定义：
捕食者-猎物动态的 Lotka-Volterra 方程
除了主要变量，还有四个参数用于描述模型中的各种生态因素：
α 表示在没有捕食者情况下猎物种群的内在增长率。β 表示捕食者对猎物的捕食率。γ 表示在没有猎物情况下捕食者种群的死亡率。δ 表示捕食者将消耗的猎物转化为新捕食者生物量的效率。
这些变量和参数共同描述了生态系统中捕食者与猎物之间的相互作用动态，并用于数学建模猎物和捕食者种群随时间的变化。这里是这个系统作为 torch.nn.Module 的实现：
class LotkaVolterra(nn.Module):
"""
The Lotka-Volterra equations are a pair of first-order, non-linear, differential equations
describing the dynamics of two species interacting in a predator-prey relationship.
"""
def __init__(self,
alpha: float = 1.5, # The alpha parameter of the Lotka-Volterra system
beta: float = 1.0, # The beta parameter of the Lotka-Volterra system
delta: float = 3.0, # The delta parameter of the Lotka-Volterra system
gamma: float = 1.0 # The gamma parameter of the Lotka-Volterra system
) -> None:
super().__init__()
self.model_params = torch.nn.Parameter(torch.tensor([alpha, beta, delta, gamma]))
def forward(self, t, state):
x = state[...,0] #variables are part of vector array u
y = state[...,1]
sol = torch.zeros_like(state)
#coefficients are part of tensor model_params
alpha, beta, delta, gamma = self.model_params
sol[...,0] = alpha*x - beta*x*y
sol[...,1] = -delta*y + gamma*x*y
return sol
def __repr__(self):
return f" alpha: {self.model_params[0].item()}, \
beta: {self.model_params[1].item()}, \
delta: {self.model_params[2].item()}, \
gamma: {self.model_params[3].item()}"
这遵循了与第一个示例相同的模式，主要区别在于我们现在有四个参数，并将它们存储为 model_params 张量。以下是捕食者-猎物方程的积分和绘图代码。
lv_model = LotkaVolterra() #use default parameters
ts = torch.linspace(0,30.0,1000)
batch_size = 30
# Create a batch of initial conditions (batch_dim, state_dim) as small perturbations around one value
initial_conditions = torch.tensor([[3,3]]) + 0.50*torch.randn((batch_size,2))
sol = odeint(lv_model, initial_conditions, ts, method='dopri5').detach().numpy()
# Check the solution
plt.plot(ts, sol[:,:,0], lw=0.5);
plt.title("Time series of the Lotka-Volterra system");
plt.xlabel("time");
plt.ylabel("x");
现在是系统的相平面图：
洛伦兹系统
我们将使用的最后一个示例是洛伦兹方程，这些方程因其美丽的混沌动态图而闻名。它们最初来源于流体动力学的简化模型，并呈现以下形式：
其中 x、y 和 z 是状态变量，σ、ρ 和 β 是系统参数。
class Lorenz(nn.Module):
"""
Define the Lorenz system as a PyTorch module.
"""
def __init__(self,
sigma: float =10.0, # The sigma parameter of the Lorenz system
rho: float=28.0, # The rho parameter of the Lorenz system
beta: float=8.0/3, # The beta parameter of the Lorenz system
):
super().__init__()
self.model_params = torch.nn.Parameter(torch.tensor([sigma, rho, beta]))
def forward(self, t, state):
x = state[...,0] #variables are part of vector array u
y = state[...,1]
z = state[...,2]
sol = torch.zeros_like(state)
sigma, rho, beta = self.model_params #coefficients are part of vector array p
sol[...,0] = sigma*(y-x)
sol[...,1] = x*(rho-z) - y
sol[...,2] = x*y - beta*z
return sol
def __repr__(self):
return f" sigma: {self.model_params[0].item()}, \
rho: {self.model_params[1].item()}, \
beta: {self.model_params[2].item()}"
这展示了如何集成这个系统并绘制结果。该系统（在这些参数值下）显示了混沌动态，因此起始条件虽然很接近，但会指数级地彼此发散。
lorenz_model = Lorenz()
ts = torch.linspace(0,50.0,3000)
batch_size = 30
# Create a batch of initial conditions (batch_dim, state_dim) as small perturbations around one value
initial_conditions = torch.tensor([[1.0,0.0,0.0]]) + 0.10*torch.randn((batch_size,3))
sol = odeint(lorenz_model, initial_conditions, ts, method='dopri5').detach().numpy()
# Check the solution
plt.plot(ts[:2000], sol[:2000,:,0], lw=0.5);
plt.title("Time series of the Lorenz system");
plt.xlabel("time");
plt.ylabel("x");
这里展示了第一组初始条件的著名蝴蝶图（相平面图）。
洛伦兹方程展示了混沌动态，并描绘出一个美丽的奇异吸引子。
数据
现在我们可以在 pytorch 中定义微分方程模型了，我们需要创建一些数据用于训练。这是事情变得真正有趣的地方，因为我们首次看到了能够利用深度学习机制来拟合参数的可能性。实际上，我们可以直接使用数据张量，但这是一种很好的组织数据的方式。如果你有一些实验数据需要使用，它也会很有用。
Torch 提供了 Dataset 类用于加载数据。要使用它，你只需创建一个子类并定义两个方法。一个是返回数据点数量的 len 函数，另一个是返回给定索引处数据点的 getitem 函数。如果你想知道这些方法是如何在 python 列表中支持 len(array) 和 array[0] 下标访问的。
其余的样板代码在父类 torch.utils.data.Dataset 中定义。我们将在定义训练循环时看到这些方法的强大功能。
class SimODEData(Dataset):
"""
A very simple dataset class for simulating ODEs
"""
def __init__(self,
ts: List[torch.Tensor], # List of time points as tensors
values: List[torch.Tensor], # List of dynamical state values (tensor) at each time point
true_model: Union[torch.nn.Module,None] = None,
) -> None:
self.ts = ts
self.values = values
self.true_model = true_model
def __len__(self) -> int:
return len(self.ts)
def __getitem__(self, index: int) -> Tuple[torch.Tensor, torch.Tensor]:
return self.ts[index], self.values[index]
接下来，让我们创建一个快速生成器函数来生成一些模拟数据以测试算法。在实际使用中，数据将从文件或数据库中加载，但在这个例子中，我们将仅生成一些数据。实际上，我建议你总是从生成的数据开始，以确保你的代码在尝试加载真实数据之前正常工作。
def create_sim_dataset(model: nn.Module, # model to simulate from
ts: torch.Tensor, # Time points to simulate for
num_samples: int = 10, # Number of samples to generate
sigma_noise: float = 0.1, # Noise level to add to the data
initial_conditions_default: torch.Tensor = torch.tensor([0.0, 0.0]), # Default initial conditions
sigma_initial_conditions: float = 0.1, # Noise level to add to the initial conditions
) -> SimODEData:
ts_list = []
states_list = []
dim = initial_conditions_default.shape[0]
for i in range(num_samples):
x0 = sigma_initial_conditions * torch.randn((1,dim)).detach() + initial_conditions_default
ys = odeint(model, x0, ts).squeeze(1).detach()
ys += sigma_noise*torch.randn_like(ys)
ys[0,:] = x0 # Set the first value to the initial condition
ts_list.append(ts)
states_list.append(ys)
return SimODEData(ts_list, states_list, true_model=model)
这只是输入一个带有一些初始状态的微分方程模型，并从中生成一些时间序列数据（并添加一些高斯噪声）。然后这些数据被传入我们的自定义数据集容器中。
训练循环
接下来我们将为 pytorch 训练循环创建一个包装函数。训练意味着我们想要更新模型参数以增加与数据的对齐度（或减少成本函数）。
深度学习中的一个技巧是不要在进行梯度步骤之前使用所有数据。这部分是因为在使用巨大数据集时，你不能将所有数据放入 GPU 的内存中，但这也可以帮助梯度下降算法避免陷入局部最小值。
训练循环的文字描述：
-
将数据集分割成小批量，这些是你整个数据集的子集。通常要随机选择这些小批量。
-
迭代小批量
-
使用当前模型参数生成预测
-
计算损失（这里我们将使用均方误差）
-
使用反向传播计算梯度。
-
使用梯度下降步骤更新参数。这里我们使用 Adam 优化器。
-
完整遍历数据集的过程称为一个周期（epoch）。
好的，这里是代码：
def train(model: torch.nn.Module, # Model to train
data: SimODEData, # Data to train on
lr: float = 1e-2, # learning rate for the Adam optimizer
epochs: int = 10, # Number of epochs to train for
batch_size: int = 5, # Batch size for training
method = 'rk4', # ODE solver to use
step_size: float = 0.10, # for fixed diffeq solver set the step size
show_every: int = 10, # How often to print the loss function message
save_plots_every: Union[int,None] = None, # save a plot of the fit, to disable make this None
model_name: str = "", #string for the model, used to reference the saved plots
*args: tuple,
**kwargs: dict
):
# Create a data loader to iterate over the data. This takes in our dataset and returns batches of data
trainloader = DataLoader(data, batch_size=batch_size, shuffle=True)
# Choose an optimizer. Adam is a good default choice as a fancy gradient descent
optimizer = torch.optim.Adam(model.parameters(), lr=lr)
# Create a loss function this computes the error between the predicted and true values
criterion = torch.nn.MSELoss()
for epoch in range(epochs):
running_loss = 0.0
for batchdata in trainloader:
optimizer.zero_grad() # reset gradients, famous gotcha in a pytorch training loop
ts, states = batchdata # unpack the data
initial_state = states[:,0,:] # grab the initial state
# Make the prediction and then flip the dimensions to be (batch, state_dim, time)
# Pytorch expects the batch dimension to be first
pred = odeint(model,
initial_state,
ts[0],
method=method,
options={'step_size': step_size}).transpose(0,1)
# Compute the loss
loss = criterion(pred, states)
# compute gradients
loss.backward()
# update parameters
optimizer.step()
running_loss += loss.item() # record loss
if epoch % show_every == 0:
print(f"Loss at {epoch}: {running_loss}")
# Use this to save plots of the fit every save_plots_every epochs
if save_plots_every is not None and epoch % save_plots_every == 0:
with torch.no_grad():
fig, ax = plot_time_series(data.true_model, model, data[0])
ax.set_title(f"Epoch: {epoch}")
fig.savefig(f"./tmp_plots/{epoch}_{model_name}_fit_plot")
plt.close()
示例
拟合 VDP 振荡器
让我们使用这个训练循环从模拟的 VDP 振荡器数据中恢复参数。
true_mu = 0.30
model_sim = VDP(mu=true_mu)
ts_data = torch.linspace(0.0,10.0,10)
data_vdp = create_sim_dataset(model_sim,
ts = ts_data,
num_samples=10,
sigma_noise=0.01)
让我们创建一个参数值错误的模型并可视化起始点。
vdp_model = VDP(mu = 0.10)
plot_time_series(model_sim,
vdp_model,
data_vdp[0],
dyn_var_idx=1,
title = "VDP Model: Before Parameter Fits");
现在，我们将使用训练循环来拟合 VDP 振荡器的参数到模拟数据中。
train(vdp_model, data_vdp, epochs=50, model_name="vdp");
print(f"After training: {vdp_model}, where the true value is {true_mu}")
print(f"Final Parameter Recovery Error: {vdp_model.mu - true_mu}")
Loss at 0: 0.1369624137878418
Loss at 10: 0.0073615191504359245
Loss at 20: 0.0009214915917254984
Loss at 30: 0.0002127257248503156
Loss at 40: 0.00019956652977271006
After training: mu: 0.3018421530723572, where the true value is 0.3
Final Parameter Recovery Error: 0.0018421411514282227
不错！让我们看看图现在是什么样子的……
图示确认我们几乎完全恢复了参数。再绘制一个图，我们绘制系统在相平面上的动态（状态变量的参数化图）。
这是该模型的训练过程的可视化：
VDP 模型的训练过程视频
洛特卡-沃尔泰拉方程
现在让我们调整方法以拟合来自洛特卡-沃尔泰拉方程的模拟数据。
model_sim_lv = LotkaVolterra(1.5,1.0,3.0,1.0)
ts_data = torch.arange(0.0, 10.0, 0.1)
data_lv = create_sim_dataset(model_sim_lv,
ts = ts_data,
num_samples=10,
sigma_noise=0.1,
initial_conditions_default=torch.tensor([2.5, 2.5]))
model_lv = LotkaVolterra(alpha=1.6, beta=1.1,delta=2.7, gamma=1.2)
plot_time_series(model_sim_lv, model_lv, data = data_lv[0], title = "Lotka Volterra: Before Fitting");
这是初始参数的拟合，然后我们将像之前一样拟合并查看结果。
train(model_lv, data_lv, epochs=60, lr=1e-2, model_name="lotkavolterra")
print(f"Fitted model: {model_lv}")
print(f"True model: {model_sim_lv}")
Loss at 0: 1.1298701763153076
Loss at 10: 0.1296287178993225
Loss at 20: 0.045993587002158165
Loss at 30: 0.02311511617153883
Loss at 40: 0.020882505923509598
Loss at 50: 0.020726025104522705
Fitted model: alpha: 1.5965800285339355, beta: 1.0465354919433594, delta: 2.817030429840088, gamma: 0.939825177192688
True model: alpha: 1.5, beta: 1.0, delta: 3.0, gamma: 1.0
首先是拟合系统的时间序列图：
现在让我们使用相平面图来可视化结果。
这是拟合过程的可视化……
洛伦兹方程
最后，让我们尝试拟合洛伦兹方程。
model_sim_lorenz = Lorenz(sigma=10.0, rho=28.0, beta=8.0/3.0)
ts_data = torch.arange(0, 10.0, 0.05)
data_lorenz = create_sim_dataset(model_sim_lorenz,
ts = ts_data,
num_samples=30,
initial_conditions_default=torch.tensor([1.0, 0.0, 0.0]),
sigma_noise=0.01,
sigma_initial_conditions=0.10)
lorenz_model = Lorenz(sigma=10.2, rho=28.2, beta=9.0/3)
fig, ax = plot_time_series(model_sim_lorenz, lorenz_model, data_lorenz[0], title="Lorenz Model: Before Fitting");
ax.set_xlim((2,15));
这是初始拟合，然后我们将调用训练循环。
train(lorenz_model,
data_lorenz,
epochs=300,
batch_size=5,
method = 'rk4',
step_size=0.05,
show_every=50,
lr = 1e-3)
这是训练结果：
Loss at 0: 114.25119400024414
Loss at 50: 4.364489555358887
Loss at 100: 2.055854558944702
Loss at 150: 1.2539702206850052
Loss at 200: 0.7839434593915939
Loss at 250: 0.5347371995449066
让我们看看拟合的模型。从完整的动态图开始。
让我们放大数据的主体部分，看看拟合效果如何。
你可以看到模型在数据范围内非常接近真实模型，并且对于 t < 16 的未见数据具有良好的泛化能力。现在是相位平面图（放大）。
plot_phase_plane(model_sim_lorenz, lorenz_model, data_lorenz[0], title = "Lorenz Model: After Fitting", time_range=(0,20.0));
你可以看到我们拟合的模型在 t ∈ [0,16] 范围内表现良好，但随后开始发散。
神经微分方程简介
这种方法在我们知道右侧方程形式的情况下效果很好，但如果我们不知道呢？我们可以使用这个方法来发现模型方程吗？
这个主题过于庞大，无法在这篇文章中完全覆盖，但将我们的微分方程模型迁移到 torch 框架的最大优势之一是我们可以将其与人工神经网络层混合搭配。
我们可以做的最简单的事情是用神经网络层替换方程右侧的 f(y,t; θ)。这类方程被称为神经微分方程，可以看作是递归神经网络的推广。
神经微分方程将方程的右侧替换为人工神经网络
作为第一个例子，让我们对我们简单的 VDP 振荡器系统进行操作。首先，我们将重新生成模拟数据，你会注意到我正在创建更长的时间序列数据和更多的样本。拟合神经微分方程需要更多的数据和计算能力，因为我们有许多需要确定的参数。
# remake the data
model_sim_vdp = VDP(mu=0.20)
ts_data = torch.linspace(0.0,30.0,100) # longer time series than the custom ode layer
data_vdp = create_sim_dataset(model_sim_vdp,
ts = ts_data,
num_samples=30, # more samples than the custom ode layer
sigma_noise=0.1,
initial_conditions_default=torch.tensor([0.50,0.10]))
现在我定义一个简单的前馈神经网络层来填充方程的右侧。
class NeuralDiffEq(nn.Module):
"""
Basic Neural ODE model
"""
def __init__(self,
dim: int = 2, # dimension of the state vector
) -> None:
super().__init__()
self.ann = nn.Sequential(torch.nn.Linear(dim, 8),
torch.nn.LeakyReLU(),
torch.nn.Linear(8, 16),
torch.nn.LeakyReLU(),
torch.nn.Linear(16, 32),
torch.nn.LeakyReLU(),
torch.nn.Linear(32, dim))
def forward(self, t, state):
return self.ann(state)
这里是拟合前系统的图：
model_vdp_nde = NeuralDiffEq(dim=2)
plot_time_series(model_sim_vdp, model_vdp_nde, data_vdp[0], title = "Neural ODE: Before Fitting");
你可以看到我们一开始离正确解很远，但我们注入到模型中的信息却少得多。让我们看看能否通过拟合模型得到更好的结果。
train(model_vdp_nde,
data_vdp,
epochs=1500,
lr=1e-3,
batch_size=5,
show_every=100,
model_name = "nde")
Loss at 0: 84.39617252349854
Loss at 100: 84.34061241149902
Loss at 200: 73.75008296966553
Loss at 300: 3.4929964542388916
Loss at 400: 1.6555403769016266
Loss at 500: 0.7814530655741692
Loss at 600: 0.41551147401332855
Loss at 700: 0.3157300055027008
Loss at 800: 0.19066352397203445
Loss at 900: 0.15869349241256714
Loss at 1000: 0.12904016114771366
Loss at 1100: 0.23840919509530067
Loss at 1200: 0.1681726910173893
Loss at 1300: 0.09865255374461412
Loss at 1400: 0.09134986530989408
通过可视化结果，我们可以看到模型能够拟合数据，甚至可以外推到未来（尽管它的效果和速度不如指定模型）。
现在是我们神经微分方程模型的相位平面图。
这些模型需要很长时间才能训练，并且需要更多数据才能获得良好的拟合。这是有道理的，因为我们同时试图学习模型和参数。
NDE 拟合过程的视频。
结论与总结
在这篇文章中，我展示了如何在 pytorch 生态系统中使用 torchdiffeq 包应用微分方程模型。这篇文章中的代码可以在github上找到，并且可以直接在google colab中打开进行实验。你也可以通过 pip 安装这篇文章中的代码：
pip install paramfittorchdemo
这篇文章是一个介绍，未来我将写更多关于以下主题的内容：
-
如何将一些动力学的机制知识与深度学习相结合。这些被称为通用微分方程，因为它们使我们能够将科学知识与深度学习结合起来。这基本上将两种方法融合在一起。
-
如何将微分方程层与其他深度学习层结合。
-
模型发现：我们能否从数据中恢复实际的模型方程？这使用SINDy等工具从数据中提取模型方程。
-
用于管理这些模型训练的 MLOps 工具。这包括MLFlow、Weights and Biases和Tensorboard等工具。
-
任何其他我从你那里听到的反馈！
如果你喜欢这篇文章，确保关注我并在lLinked-in上与我联系。除非另有说明，否则所有图片均由作者提供。
扩散概率模型与文本到图像生成
原文：
towardsdatascience.com/diffusion-probabilistic-models-and-text-to-image-generation-9f441d0bc786
你能想象的任何东西的照片级生成
·发布于 Towards Data Science ·阅读时间 4 分钟·2023 年 3 月 29 日
–
图 1. 文本到图像生成。图片由作者制作。
如果你是最新计算机视觉论文的忠实追随者，你一定会对生成网络在图像生成中的惊人效果感到惊讶。许多之前的文献都是基于开创性的生成对抗网络（GAN）思想，但最近的论文情况已不再如此。事实上，如果你仔细查看最新的论文，如 ImageN 和 Staple Diffusion，你会不断看到一个陌生的术语：扩散概率模型。
本文深入探讨了新兴趋势模型的基本知识，简要概述了如何学习，以及随之而来的令人兴奋的应用。
从逐渐添加高斯噪声开始……
图 2. 扩散去噪概率模型概述。图片来源：arxiv.org/abs/2006.11239。
考虑一张添加了少量高斯噪声的图像。图像可能变得有些嘈杂，但原始内容仍然很可能被识别出来。现在重复这一步骤；最终，图像将变成几乎纯粹的高斯噪声。这被称为扩散概率模型的前向过程。
目标很简单：通过利用前向过程是马尔可夫链的事实（当前时间框架的过程与之前的时间框架独立），我们实际上可以学习一个逆过程，在当前帧上稍微去噪图像。
给定一个经过适当学习的逆过程和一个随机高斯噪声，我们现在可以反复应用噪声，并最终获得一个与原始数据分布非常相似的图像——因此这是一个生成模型。
扩散模型的一个优点是，训练可以通过仅选择中间的随机时间戳进行优化（而不是必须完全重建图像）。与 GAN 相比，训练本身更加稳定，因为小的超参数差异可能会导致模型崩溃。
请注意，这只是对去噪扩散概率模型的一个非常高层次的概述。有关数学细节，请参考 这里 和 这里*。
文字到图像生成的惊人结果
图 3. ImageN 生成的结果。文本提示位于图像下方。图像来源于：arxiv.org/abs/2205.11487。
图像生成中的去噪扩散模型的理念首次发布于 2020 年，但直到最近的 Google 论文ImageN才真正引爆了这个领域。
类似于 GAN，扩散模型也可以根据图像和文本等提示进行条件化。Google Research Brain 团队建议，冻结的大型语言模型实际上是提供文本条件以生成真实感图像的绝佳编码器。
从 2D 到 3D 的转变……
图 4. DreamFusion 流程概述。图像来源于：arxiv.org/abs/2209.14988。
与众多计算机视觉趋势一样，在二维领域的出色表现激发了向 3D 扩展的雄心；扩散模型也不例外。最近，Poole 等人提出了基于 ImageN 和 NeRF 坚实基础的DreamFusion文本到 3D 模型。
有关 NeRF 的简要概述，请参阅 这里*。
图 4 展示了 DreamFusion 的流程图。该流程从一个随机初始化的 NeRF 开始。基于生成的密度、反射率和法线（在给定光源下），网络输出着色，并随后根据特定摄像机角度的 NeRF 形成颜色。渲染的图像与高斯噪声结合，目标是利用冻结的 ImageN 模型重建图像，并随后更新 NeRF 模型。
图 5. DreamFusion 的结果。图像来源于：arxiv.org/abs/2209.14988。
一些令人惊叹的 3D 结果展示在图 5 所示的画廊中。通过简单图像完全展现了对象的一致颜色和形状。
最近的工作如Magic3D进一步改进了流程，使重建过程更快且更细致。
结束备注
这就是对图像生成扩散模型进展的概述。当简单的词汇转变为生动的图像时，大家更容易想象并描绘他们最疯狂的想法。
“写作是声音的绘画” — 伏尔泰
感谢你能看到这里 🙏*!* 我定期撰写有关计算机视觉/深度学习的不同领域的文章，如果你有兴趣了解更多， 请加入并订阅 ！
使用 Python 和 MySQL 进行数字营销分析
原文：
towardsdatascience.com/digital-marketing-analysis-simultaneously-with-python-and-mysql-ee00e05a3813
一个数字营销分析练习，展示了 SQL 和 Python 环境中的逐步解释代码
·发布于Towards Data Science ·阅读时间 15 分钟·2023 年 3 月 19 日
–
照片由Zdeněk Macháček拍摄，发布在Unsplash上
介绍
在这段简短的旅程中，我们将探索一个简单的数据集，包含基本的网页营销指标，如‘用户’，‘会话’和‘跳出率’，为期五个月。
这种设置的目的是获取一些基本但有用的知识，以回答一些必要的操作营销问题，而不是专注于理解网站性能。
我们将关注两种强大且常用的数字工具，探索两种方法，以便最终得出相同的结果。
一方面，我们将探讨MySQL Workbench的语法，涉及一些多样化的查询；另一方面，对于每个问题，我们将使用Python的图形和可视化资源进行对比。这两个环境将分别标记为# MySQL 和# Python。每个问题将附有注释和两种代码的解释，以便更深入理解。
# MySQL
-- displaying dataset (case_sql.csv)
SELECT * FROM case_sql;
图片由作者提供。
你可以在这里下载 SQL 数据集。
# Python
# import python libraries
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.style as style
%matplotlib inline
color = sns.color_palette()
from pandas.plotting import table
from datetime import datetime
# load the data set
df = pd.read_csv("case.csv", sep=";")
# number of observations and columns
df.shape
(31507, 7)
# display rows sample
df.sample(15)
图片由作者提供。
# SHAPE
# Create a function that prints the shape of the dataframe and some other basic info
# --> number of observations, features, duplicates, missing values (True, False) datatypes and its proportion.
def shape_df(df):
print(f"Number of observations: {df.shape[0]}")
print(f"Number of variables: {df.shape[1]}")
print(f"Number of duplicates: {df.duplicated().sum()}")
print(f"Are there any missing values? {df.isnull().values.any()}\n-----")
print(f"{df.dtypes.sort_values(ascending=True)}\n-----")
print(f"Datatypes' proportion:\n{df.dtypes.value_counts(ascending=True)}")
# calling the function
shape_df(df)
Number of observations: 31507
Number of variables: 7
Number of duplicates: 4083
Are there any missing values? False
-----
date int64
users int64
sessions int64
bounces int64
brand object
device_category object
default_channel_grouping object
dtype: object
------
Datatypes proportion:
object 3
int64 4
dtype: int64
# lowering columns' capital letters for easy typing
df.columns = map(str.lower, df.columns)
# remove spaces in columns name
df.columns = df.columns.str.replace(' ','_')
将日期转换为 datetime 类型
# make string version of original column 'date', call it 'date_'
df['date_'] = df['date'].astype(str)
# create the new columns using string indexing
df['year'] = df['date_'].str[0:4]
df['month'] = df['date_'].str[4:6]
df['day'] = df['date_'].str[6:]
# concatenate 'year', 'month' and 'day'
df["date"] = df["year"] + "-" + df["month"] + "-" + df["day"]
# convert to datetime
df["date"] = pd.to_datetime(df["date"], format="%Y-%m-%d")
# extract 'year', 'month' and 'weekday'
df["year"] = df["date"].dt.year
df["month"] = df["date"].dt.month
df["weekday"] = df["date"].dt.dayofweek.map({0 : "Mon", 1 : "Tue", 2 : "Wed", 3: "Thu", 4 : "Fri", 5 : "Sat", 6 : "Sun"})
# select columns to perform exploratory data analysis
cols = "date year month weekday brand device_category default_channel_grouping users sessions bounces".split()
df = df[cols].copy()
# display final dataset
df.head(10)
图片由作者提供。
你可以在这里下载 Python 数据集。
1. 查看按设备类型划分的用户分布
# MySQL
-- Device distribution
SELECT
device_category, -- select the device_category column
ROUND(COUNT(users) / (SELECT
COUNT(users)
FROM
case_sql) * 100,
1) AS percent -- calculate the percentage of users in each category
FROM
case_sql -- select data from the case_sql table
GROUP BY 1 -- group the result by device_category
ORDER BY 1; -- order the result by device_category in ascending order
图片来源于作者。
我们可以看到，通过设备分布，手机和桌面并列为最常见的访问类型。
# Python
# device distribution
# counts the number of occurrences of each unique device category in the device_category column of the DataFrame df, including missing values (if any).
df.device_category.value_counts(dropna=False).plot(kind='pie', figsize=(8,4),
explode = (0.02, 0.02, 0.02),
autopct='%1.1f%%',
startangle=150);
# Params
plt.xticks(rotation=0, horizontalalignment="center")
plt.title("Device distribution", fontsize=10, loc="right");
图片来源于作者。
2. 从品牌的角度来看，上述情况如何？
# MySQL
-- Brand distribution
SELECT
brand, -- select the brand column
COUNT(users) AS users, -- count the number of users for each unique brand and alias the result as "users"
ROUND(COUNT(users) / (SELECT
COUNT(users)
FROM
case_sql) * 100,
2) AS percent -- calculate the percentage of users for each brand out of the total number of users in the case_sql table and alias the result as "percent"
FROM
case_sql -- select data from the case_sql table
GROUP BY 1; -- group the result by the first column (brand)
图片来源于作者。
在品牌层面上，品牌 2 的访问量占比最高，为 56.28%，而品牌 1 的总访问量为 43.72%。
# Python
# Brand distribution
absolut = df["brand"].value_counts().to_frame()
# Pie chart
absolut.plot(kind='pie', subplots=True, autopct='%1.2f%%',
explode= (0.05, 0.05), startangle=20,
legend=False, fontsize=12, figsize=(8,4))
# Params
plt.xticks(rotation=0, horizontalalignment="center")
plt.title("Brand's distribution", fontsize=10, loc="right");
display(absolut) # Table
图片来源于作者。
3. 我们看到最多用户访问品牌 1 网站的那一天是星期几？
# MySQL
SELECT
date, -- select the date column
DAYNAME(date) AS day_name, -- calculate the day name corresponding to each date
SUM(users) AS users -- sum the number of users for each unique date where the brand is 'Brand 1'
FROM
case_sql -- select data from the case_sql table
WHERE
brand = 'Brand 1' -- filter rows where the brand is 'Brand 1'
GROUP BY 1 -- group the result by date
ORDER BY 3 DESC -- order the result by users in descending order
LIMIT 1; -- select only the first row of the result (the row with the highest number of users)
图片来源于作者。
在 2019 年 9 月至 2020 年 1 月之间的 298,412 名用户中，访问品牌 1网站的用户最多的一天是2019 年 11 月 22 日，共计 885 次访问，那天是周五。
# Python
# filter users that arrived at 'Brand 1' only, assign it 'brand_1'
brand_1 = df[df["brand"] == "Brand 1"].copy()
''' sum total users that came from all "channelgrouping" for the same date,
assign it 'brandgroup' no matter the type of device '''
brandgroup = brand_1.groupby(["date","weekday"])[["default_channel_grouping","users"]].sum()
# filter the date by maximum users, assign it 'users'
users = brandgroup[brandgroup["users"] == brandgroup.users.max()].copy()
# reseat index
users.reset_index(["date"], inplace=True)
users.reset_index(["weekday"], inplace=True)
# results
print(f"""Date: {users.date} \n\nTotal users: {users.users} \n\nDay of week: {users.weekday}""")
Date: 0 2019-11-22
Name: date, dtype: datetime64[ns]
Total users: 0 885
Name: users, dtype: int64
Day of week: 0 Fri
Name: weekday, dtype: object
# calling the variable
users
图片来源于作者。
3.1 在那一天，有多少用户访问了品牌 2？
# MySQL
SELECT
DATE(date) AS date, -- Select the date from the 'date' column and convert it to a date data type
DAYNAME(date) AS dayofweek, -- Select the day of the week from the 'date' column
SUM(CASE
WHEN brand = 'Brand 1' THEN users -- Sum the 'users' column for Brand 1
ELSE NULL
END) AS b1_users,
SUM(CASE
WHEN brand = 'Brand 2' THEN users -- Sum the 'users' column for Brand 2
ELSE NULL
END) AS b2_users
FROM
case_sql -- From the 'case_sql' table
GROUP BY 1, 2 -- Group the results by the first and second columns (date and dayofweek)
ORDER BY 3 DESC -- Order the results by b1_users in descending order
LIMIT 1; -- Limit the results to only the highest total number of Brand 1 users
图片来源于作者。
实际上，两个品牌在同一天都达到了最高访问量。
# Python
# filter users that arrived at 'Brand 2', assign it 'brand_2'
brand_2 = df[df["brand"] == "Brand 2"].copy()
# rename the 'users' column from previous (above) Python code
brandgroup.rename(columns = {'users':'brand1_users'}, inplace = True)
# include a new column with the filtered users from 'Brand_2'
brandgroup["brand2_users"] = brand_2.groupby(["date","weekday"])[["default_channel_grouping","users"]].sum()
# filter the new column (brand2_users) by maximum users
users2 = brandgroup[brandgroup["brand2_users"] == brandgroup.brand2_users.max()].copy()
图片来源于作者。
4. 在所有渠道分组中，哪个渠道贡献了最高数量的用户？
# MySQL
SELECT
default_channel_grouping AS channels,
SUM(users) AS total_users,
ROUND(SUM(users) / (SELECT
SUM(users)
FROM
case_sql) * 100,
1) AS percent -- calculate the percentage of users for each channel
FROM
case_sql
GROUP BY 1
ORDER BY 2 DESC;
图片来源于作者。
自然搜索仍然是生成最多用户的渠道（近 141,000 名用户），占据了两个网站总访问量的近一半，其次是付费搜索和直接访问。展示广告排名第 4，社交媒体排名第 6，共贡献了 6,722 位用户。
# Python
# sum users by all channel groups and plot bar chart
ax = df.groupby("default_channel_grouping")["users"].sum().sort_values(ascending=True)\
.plot(kind="bar", figsize=(9,6), fontsize=12, linewidth=2,
color=sns.color_palette("rocket"), grid=False, table=False)
# show data labels
for p in ax.patches:
ax.annotate("%.0f" % p.get_height(), (p.get_x() + p.get_width() / 2., p.get_height()),
ha='center', va='center', xytext=(0, 7), textcoords='offset points')
# params
plt.xlabel("Channel groups", fontsize=10)
plt.xticks(rotation=90, horizontalalignment="center")
plt.ylabel("Absolute values", fontsize=10)
plt.title("Best channel group (highest number of users)", fontsize=10, loc="right");
图片来源于作者。
4.1 所有渠道分组中，哪个品牌贡献的用户数量最高？
# MySQL
SELECT
default_channel_grouping AS channels,
SUM(CASE -- sum users by brand and map to new columns
WHEN brand = 'brand 1' THEN users -- if brand = 'brand 1', sum users and store in 'Brand_1' column
ELSE NULL -- if not 'brand 1', set value to null
END) AS Brand_1, -- create column for Brand 1 users
SUM(CASE
WHEN brand = 'brand 2' THEN users
ELSE NULL
END) AS Brand_2
FROM
case_sql
GROUP BY 1 -- group by channel
ORDER BY 3 DESC; -- order by Brand 2 users in descending order
图片来源于作者。
# Python
# create pivot_table
# sum all users for each brand by channels
type_pivot = df.pivot_table(
columns="brand",
index="default_channel_grouping",
values="users", aggfunc=sum)
display(type_pivot)
#Display pivot_table with a bar chart
type_pivot.sort_values(by=["Brand 2"], ascending=True).plot(kind="bar", figsize=(12,8) ,fontsize = 15)
plt.xlabel("Channel groups", fontsize=10)
plt.xticks(rotation=90, horizontalalignment="center")
plt.ylabel("Absolute values", fontsize=10)
plt.title("Channel groups by brand (highest number of users)", fontsize=10, loc="right");
图片来源于作者。
图片来源于作者。
自然搜索为品牌 2 带来了 105,062 位用户，为品牌 1 带来了 35,911 位用户。除‘其他’之外，品牌 1 在这一方面表现更优，品牌 2 在所有渠道中为网站带来的用户数量最高。
5. 在所有渠道中，哪个品牌在 2019 年贡献了至少 5% 的付费会话？
# MySQL
SELECT
brand,
default_channel_grouping AS channels,
ROUND(SUM(sessions) / (SELECT
SUM(sessions)
FROM
case_sql) * 100,
1) AS percent
FROM
case_sql
WHERE
default_channel_grouping IN ('Paid Search' , 'Paid Social', 'Display', 'Other Advertising') -- include only rows with these values
AND date < '2020-01-01' -- only date before '2020-01-01' will be included.
GROUP BY 1 , 2
HAVING percent > 5 -- filters the groups to only include values greater than 5%.
ORDER BY 1 , 3 DESC
图片来源于作者
# Python
# groupby dataframe by selected cols
df = df.groupby(["date","brand","default_channel_grouping"])["sessions"].sum().to_frame().copy()
# calculate percentages (new column)
df["percent"] = (df.apply(lambda x: x/x.sum())*100).round(2)
# reset index
df = df.reset_index().copy()
# display a 5 rows sample
df.sample(5)
图片来源于作者。
# filter paid channels using lambda function
paid = df.apply(lambda row: row[df['default_channel_grouping'].isin(['Display','Paid Search','Paid Social','Other Advertising'])])
# filter year 2019
paid = paid[paid['date'] < '2020-01-01']
# groupby channels by brand
paid = paid.groupby(["brand","default_channel_grouping"])[["sessions","percent"]].sum()
# filter sessions higher than 5%
paid[paid["percent"] >5]
图片来源于作者。
6. 两个品牌按设备类型分别收到多少次访问？
# MySQL
SELECT
brand,
SUM(CASE
WHEN device_category = 'Desktop' THEN users
ELSE NULL
END) AS desktop,
SUM(CASE
WHEN device_category = 'Mobile' THEN users
ELSE NULL
END) AS mobile,
SUM(CASE
WHEN device_category = 'Tablet' THEN users
ELSE NULL
END) AS tablet
FROM
case_sql
GROUP BY 1
ORDER BY 1;
图片来源于作者。
# Python
# pivot_table
type_pivot = df.pivot_table(
columns="device_category",
index="brand",
values="users", aggfunc=sum)
display(type_pivot)
# display pivot_table (chart)
ax = type_pivot.sort_values(by=["brand"], ascending=False).plot(kind="bar", figsize=(12,8) ,fontsize = 15);
# adding data labels
for p in ax.patches:
ax.annotate("%.0f" % p.get_height(), (p.get_x() + p.get_width() / 2., p.get_height()), ha='center', va='center', xytext=(0, 7), textcoords='offset points')
plt.xlabel("Brands", fontsize=10)
plt.xticks(rotation=0, horizontalalignment="center")
plt.ylabel("Absolute values", fontsize=10)
plt.title("Brand by type of device", fontsize=10, loc="right")
plt.legend(["Desktop","Mobile","Tablet"]);
图片来源于作者。
移动设备 是品牌 2 的首选设备类型，而桌面设备 是品牌 1 使用最多的设备。
6.1 用户平均设备使用类型按渠道分布如何？
# MySQL
SELECT
default_channel_grouping,
AVG(CASE
WHEN device_category = 'Desktop' THEN users
ELSE NULL
END) AS desktop,
AVG(CASE
WHEN device_category = 'Mobile' THEN users
ELSE NULL
END) AS mobile,
AVG(CASE
WHEN device_category = 'Tablet' THEN users
ELSE NULL
END) AS tablet
FROM
case_sql
GROUP BY 1
ORDER BY 1;
图片来自作者。
# Python
# pivot_table
type_pivot = df.pivot_table(
columns="device_category",
index="default_channel_grouping",
values="users", aggfunc=np.mean)
display(type_pivot)
图片来自作者。
# display pivot_table
type_pivot.sort_values(by=["default_channel_grouping"], ascending=False).plot(kind="bar", figsize=(12,8) ,fontsize = 15);
plt.xlabel("Date", fontsize=10)
plt.xticks(rotation=90, horizontalalignment="center")
plt.ylabel("Absolute values", fontsize=10)
plt.title("Average use of device types by channel grouping", fontsize=10, loc="right")
plt.legend(["Desktop","Mobile","Tablet"]);
图片来自作者。
平均而言，桌面设备在推荐、直接和其他渠道中使用更频繁。至于其他渠道，垂直导向的内容应该始终考虑在内。
7. 如何评估渠道分组的跳出率？
跳出率的计算方法是将总跳出次数除以总会话次数。
# MySQL
SELECT
default_channel_grouping,
SUM(sessions) AS sessions,
SUM(bounces) AS bounces,
ROUND(SUM(bounces) / SUM(sessions) * 100, 2) AS bounces_r
FROM
case_sql
GROUP BY 1
ORDER BY 4 DESC;
图片来自作者。
平均跳出率：54.93% (avg_bounces_r)
SELECT
SUM(sessions) AS sessions,
SUM(bounces) AS bounces,
ROUND(SUM(bounces) / SUM(sessions) * 100, 2) AS bounces_r,
AVG(ROUND(bounces/sessions*100, 2)) AS avg_bounces_r
FROM
case_sql;
# Python
# group individual channels by sum of users
dfbounce = df.groupby("default_channel_grouping")["users"].sum().to_frame()
# group individual channels by sum of sessions
dfbounce["sessions"] = df.groupby("default_channel_grouping")["sessions"].sum()
# group individual channels by sum of bounces
dfbounce["bounces"] = df.groupby("default_channel_grouping")["bounces"].sum()
# calculus of bounce rate for each individual channel
dfbounce["bounces_r"] = dfbounce.apply(lambda x: 0.0 if x["sessions"] == 0.0 else (x["bounces"] / x["sessions"])*100, axis=1).round(2)
dff = dfbounce.copy()
dfbounce.drop(["users"],axis=1,inplace=True)
# sort values by rate
dfbounce.sort_values(by="bounces_r", ascending=False)
图片来自作者。
# display bar chart with the bounce rate for each channel
ax = dfbounce.groupby("default_channel_grouping")["bounces_r"].sum().sort_values(ascending=True)\
.plot(kind="bar", figsize=(9,6), fontsize=12, linewidth=2, color=sns.color_palette("rocket"), grid=False, table=False)
for p in ax.patches:
ax.annotate("%.2f" % p.get_height(), (p.get_x() + p.get_width() / 2., p.get_height()),
ha='center', va='center', xytext=(0, 7), textcoords='offset points')
plt.axhline(dfbounce.groupby("default_channel_grouping")["bounces_r"].mean().mean(), linewidth=1, color ="r")
plt.xlabel("channel groups", fontsize=10)
plt.xticks(rotation=90, horizontalalignment="center")
plt.ylabel("Absolute values", fontsize=10)
plt.title("Bounce rate by channelGrouping", fontsize=10, loc="right");
图片来自作者。
正如预期，显示渠道的跳出率最高，其次是直接渠道和付费社交。自然搜索的跳出率与平均水平持平。低于最低跳出率阈值的有推荐、原生广告和其他广告。
7.1 网站上的跳出率是否随着时间的推移有所改善或恶化？
# MySQL
SELECT
YEAR(date) AS year, -- extract year
MONTH(date) AS month, -- extract month
DATE_FORMAT(date, '%b') AS month_, -- format the date column to display month name
ROUND(SUM(bounces) / SUM(sessions) * 100, 2) AS bounces_r -- calculate bounce rate
case_sql
GROUP BY 1 , 2 , 3
ORDER BY 1 , 2 , 3;
图片来自作者。
# Python
df_date = df.groupby("date")[['sessions','bounces']].sum()
''' create function to assess the bounce rate, assign it as 'bounce_r'
Return 0 if session's value is 0, else divide the bounces by sessions
for each date and multiply it by 100 to get the percentage '''
def div(bounces, sessions):
return lambda row: 0.0 if row[sessions] == 0.0 else float((row[bounces]/(row[sessions])))*100
# create column 'bounce_r' with the function results
df_date["bounce_r"] = (df_date.apply(div('bounces', 'sessions'), axis=1)).round(1)
# drop unnecessary columns
df_date.drop(["sessions","bounces"], axis=1, inplace=True)
# sum all bounces over time and plot chart
ax = df_date.plot(kind="line", figsize=(14,6), fontsize=12, linewidth=2)
plt.xlabel("Date", fontsize=10)
plt.xticks(rotation=90, horizontalalignment="center")
plt.ylabel("Rate", fontsize=10)
plt.title("Evolution of the bounce rate over time", fontsize=10, loc="right");
图片来自作者。
# Smoothing the line with a step of 15 days interval
resampled = df_date["bounce_r"].resample("m").mean()
plt.figure(figsize = (12,6))
ax = sns.lineplot(data = resampled)
plt.title("Evolution of the bounce rate over time (smooth)", fontsize=10, loc="right")
plt.xlabel("Date", fontsize=10)
plt.xticks(rotation=0, horizontalalignment="center")
plt.ylabel("Rate", fontsize=10);
图片来自作者。
网站上的跳出率正在改善。
7.2 跳出率按渠道和品牌的细分
# Python
# filter by brand
b1 = df[df["brand"] == "Brand 1"]
b2 = df[df["brand"] == "Brand 2"]
# ** brand 1 **
# group individual channels by sum of sessions for brand 1
dfbrand = b1.groupby("default_channel_grouping")["sessions"].sum().to_frame()
dfbrand.rename(columns={"sessions":"sessions1"}, inplace=True)
# group individual channels by sum of bounces for brand 1
dfbrand["bounces1"] = b1.groupby("default_channel_grouping")["bounces"].sum()
# calculus of bounce rate for each individual channel for brand 1
dfbrand["1bounces_r"] = dfbrand.apply(lambda x: 0.0 if x["sessions1"] == 0.0 else (x["bounces1"] / x["sessions1"]*100), axis=1).round(2)
# ** brand 2 **
# group individual channels by sum of bounces for brand 2
dfbrand["sessions2"] = b2.groupby("default_channel_grouping")["sessions"].sum()
# group individual channels by sum of bounces for brand 2
dfbrand["bounces2"] = b2.groupby("default_channel_grouping")["bounces"].sum()
# calculus of bounce rate for each individual channel for brand 2
dfbrand["2bounces_r"] = dfbrand.apply(lambda x: 0.0 if x["sessions2"] == 0.0 else (x["bounces2"] / x["sessions2"]*100), axis=1).round(2)
# sort values by rate
dfbrand.sort_values(by="1bounces_r", ascending=False)
图片来自作者。
# clean dataframe
dfchannels = dfbrand.copy()
dfbrand_chart = dfbrand.copy()
dfbrand_chart.drop(["sessions1","sessions2","bounces1","bounces2"], axis=1, inplace=True)
# display bar chart with the average bounce rate for each channel
ax = dfbrand_chart.plot(kind="bar", figsize=(13,6), fontsize=12, linewidth=2, color=sns.color_palette("BrBG"), grid=False, table=False)
for p in ax.patches:
ax.annotate("%.1f" % p.get_height(), (p.get_x() + p.get_width() / 2., p.get_height()), ha='center', va='center', xytext=(0, 7), textcoords='offset points')
plt.xlabel("channel groups", fontsize=10)
plt.xticks(rotation=90, horizontalalignment="center")
plt.ylabel("Absolute values", fontsize=10)
plt.title("Bounce rate by channelGrouping and by brand", fontsize=10, loc="right")
plt.legend(["Brand 1","Brand 2"]);
图片来自作者。
品牌 1 的显示渠道跳出率低于预期。通常，这些渠道会导致较高的跳出率。需要详细了解内容策略，并尝试将其调整到品牌 2。
品牌 2 在推荐渠道中的跳出率高于可接受范围。
8. 入站媒体与付费媒体的比例
# MySQL
SELECT
brand,
CASE
WHEN
default_channel_grouping IN ('Paid Search',
'Paid Social',
'Display',
'Other Advertising')
THEN
'Paid'
WHEN
default_channel_grouping IN ('Direct',
'Native',
'Organic Search',
'Referral',
'Social',
'Email',
'(Other)')
THEN
'Organic'
ELSE NULL
END AS media,
ROUND(SUM(bounces) / SUM(sessions) * 100, 2) AS bounce_r
FROM
case_sql
GROUP BY brand , media
ORDER BY 1;
图片来自作者。
# Python
# create dictionary
media_dict =
{
'Display': 'paid',
'Paid Search': 'paid',
'Paid Social': 'paid',
'Other Advertising': 'paid',
'Direct': 'organic',
'Native': 'organic',
'Organic Search': 'organic',
'Referral': 'organic',
'Social': 'organic',
'Email': 'organic',
'(Other)': 'organic'
}
# mapping the dict into a new column
df['media'] = df['default_channel_grouping'].map(media_dict)
# define cols position in dataframe
cols = ['brand','media','sessions','bounces']
# reindex columns order
df = df.reindex(columns = cols)
# groupby dataframe by selected cols
df = df.groupby(["brand","media"])[["sessions","bounces"]].sum()
# bounce rate by channel
df["bounces_r"] = df.apply(lambda x: 0.0 if x["sessions"] == 0.0 else (x["bounces"] / x["sessions"])*100, axis=1).round(2)
图片来自作者。
结论
正如承诺的，我们通过逐步的方法进行了简单的数字营销分析，使用了 MySQL Workbench 和 Python。
两种工具各有其特点和要求，但推理过程相对类似，抛开它们的图形能力和局限性。
欢迎下载数据集并通过实践本文中涉及的一些技术细节来探索，实施新代码并提出进一步的分析问题。
探索你可能也喜欢的其他项目：
在 MySQL Workbench 中审查一些应用于营销操作和常见问题的主要 SQL 查询的方法……
towardsdatascience.com ## 15 个关于移动营销活动的商业问题：ROAS（广告支出回报）
一项探索性营销数据分析，用于监控和评估移动营销活动的表现（EDA）
towardsdatascience.com ## 机器学习：预测银行贷款违约
一种数据科学方法，用于预测和理解申请人的档案，以最小化未来贷款违约的风险。
towardsdatascience.com
如何与我联系：
✅ 感谢阅读！
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
