使用量子退火进行 scikit-learn 特征选择
对于具有大量特征的数据集,使用量子处理进行 scikit-learn 模型的特征选择
·
关注 发表在 Towards Data Science · 11 分钟阅读 · 2023 年 4 月 10 日
–
特征选择是机器学习中的一个广泛话题。正确实施时,它可以帮助减少过拟合、提高可解释性、降低计算负担等。用于特征选择的技术有很多种,它们的共同点在于它们会查看特征集,并尝试将那些能够带来良好结果(准确的预测、可解释的模型等)的特征与那些不利于此目标的特征分开。
特别困难的是特征数量非常大的情况。对所有特征组合进行穷举探索通常计算开销很大。像 R 中的regsubsets()
函数那样的逐步搜索算法可能会通过添加/删除特征并比较结果来尝试推断出有前途的特征组合。但最终,当特征数量很大时,搜索成功的程度与寻找最佳特征组合所花费的努力之间仍然存在权衡。
一般来说,当问题的最佳解涉及到大量组合的搜索时,量子退火可能值得研究。我将展示一个使用 D-Wave 最近发布的 scikit-learn 插件进行特征选择的例子,该数据集具有数百个特征。
D-Wave 和 scikit-learn
请记住,这并不是通用的门模型量子计算。这是一种算法,本质上类似于模拟退火,其中有一个目标函数,并且使用类似于模拟退火的方法来找到一个最小化目标的值组合。不同的是,这里的退火并非模拟的,而是通过编程使实际系统的物理能量与目标函数相匹配。系统的能量会降低直到达到最低点(退火),然后解就是系统的状态,这个状态会被读取并返回给用户。
D-Wave 构建了量子退火器,可以有效地解决许多组合复杂度很高的问题。只要你能将问题简化为二次二进制模型(BQM),或者带约束的 BQM(CQM),或者上述模型的某些离散推广(DQM),问题就可以提交给量子解算器。更多细节请参见文档。
理论上,你可以将特征选择算法表述为 BQM,其中特征的存在是值为 1 的二进制变量,而特征的缺失是值为 0 的变量,但这需要一些工作。D-Wave 提供了一个可以直接插入到 scikit-learn 管道中的scikit-learn 插件,简化了这一过程。
本文将首先展示显式地表述 BQM 的整个过程,然后是 CQM,将模型发送到量子解算器,然后解析结果。这是解决量子退火器上优化过程的一般方法,并且可以用于许多问题。
但如果你只是想选择数据集中最好的特征,一个简单的SelectFromQuadraticModel()
方法调用就足够了。这将整个算法压缩成一行代码。我们将在最后展示这一点。
问题与方法
考虑从 OpenML 下载的场景识别数据集,最初由Boutell, M., Luo, J., Shen, X., Brown, C. (2004)创建,并由其作者在 CC BY 许可证下分发。这是一个二分类数据集,具有一个依赖变量(Urban),并有近 300 个特征。响应变量的值(二元:城市或非城市)需要基于特征的组合进行预测。一个简单的分类模型如RandomForestClassifier()
应该表现良好,前提是特征选择得当。
>>> dataset.get_data()[0]
attr1 attr2 attr3 attr4 attr5 attr6 attr7 attr8 ... attr293 attr294 Beach Sunset FallFoliage Field Mountain Urban
0 0.646467 0.666435 0.685047 0.699053 0.652746 0.407864 0.150309 0.535193 ... 0.014025 0.029709 1 0 0 0 1 0
1 0.770156 0.767255 0.761053 0.745630 0.742231 0.688086 0.708416 0.757351 ... 0.082672 0.036320 1 0 0 0 0 1
2 0.793984 0.772096 0.761820 0.762213 0.740569 0.734361 0.722677 0.849128 ... 0.112506 0.083924 1 0 0 0 0 0
3 0.938563 0.949260 0.955621 0.966743 0.968649 0.869619 0.696925 0.953460 ... 0.049780 0.090959 1 0 0 0 0 0
4 0.512130 0.524684 0.520020 0.504467 0.471209 0.417654 0.364292 0.562266 ... 0.164270 0.184290 1 0 0 0 0 0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2402 0.875782 0.901653 0.926227 0.721366 0.795826 0.867642 0.794125 0.899067 ... 0.254413 0.134350 0 0 0 0 0 1
2403 0.657706 0.669877 0.692338 0.713920 0.727374 0.750354 0.684372 0.718770 ... 0.048747 0.041638 0 0 0 0 0 1
2404 0.952281 0.944987 0.905556 0.836604 0.875916 0.957034 0.953938 0.967956 ... 0.017547 0.019734 0 0 0 0 0 1
2405 0.883990 0.899004 0.901019 0.904298 0.846402 0.858145 0.851362 0.852472 ... 0.226332 0.223070 0 0 0 0 0 1
2406 0.974915 0.866425 0.818144 0.936140 0.938583 0.935087 0.930597 1.000000 ... 0.025059 0.004033 0 0 0 0 0 1
[2407 rows x 300 columns]
一种选择特征的方法,描述在Milne, Rounds, 和 Goddard (2018)的论文中,适合量子退火机,并可以总结如下:
将数据集拆分为特征矩阵和响应向量——即 scikit-learn 中熟悉的(X, y)元组。使用这些,构造一个相关矩阵 C,其中 Cii 表示特征与响应的相关性,Cij 是特征之间的相互相关性。对于与响应高度相关的特征,我们希望尽可能多地捕获。对于彼此相关的特征,我们希望尽可能少地捕获,但不影响与响应强相关的特征。显然,我们希望在所有这些标准之间取得平衡。
设 Xi 为二元变量,指示特征 i 是否被选择到最终数据集中。如果选择了特征 i,则 Xi = 1,否则 Xi = 0。问题变成了等同于寻找极值:
目标函数
第一个项的总和表示来自特征的个别贡献——我们称之为线性项。第二个项的总和可以说包含了二次交互项。alpha 是一个偏置系数,它控制目标函数中我们允许的特征之间的交互量;其值需要探索以找到最佳结果。找到最小化目标函数的 Xi 集合等同于特征选择。
目标函数实际上是一个二次二元模型——BQM。它是二元的,因为 Xi 可以是 0 或 1。它是二次的,因为最高阶的项是二次交互项。这可以很容易地在量子退火机上解决。我们需要应用的唯一约束是,等于 1 的变量 Xi 的数量不能超过我们愿意接受的特征总数。
用困难的方法进行特征选择
让我们来解决这个问题。下面的代码块导入了我们需要的所有库,下载了数据集,并实例化了一个分类模型。
import numpy as np
import openml
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
import plotly.express as px
from plotly.subplots import make_subplots
import dimod
from dwave.system import LeapHybridCQMSampler
dataset = openml.datasets.get_dataset(312)
X, y, categorical_indicator, attribute_names = dataset.get_data(
target=dataset.default_target_attribute, dataset_format="dataframe"
)
X = X.astype(float)
y = y.values.astype(int)
clf = RandomForestClassifier()
我们关注两件事:
-
每个特征的“相关性”,即它与响应的相关强度
-
特征的“冗余性”,即相关矩阵中的非对角项(二次交互项的权重)
让我们构建几个函数来可视化相关性和冗余性,将它们应用于没有特征选择的整个数据集,然后在整个数据集上进行交叉验证(5 折)模型的性能。
def evaluate_model(m, X, y, indices=None):
if not indices is None:
X_filtered = X.iloc[:, indices]
else:
X_filtered = X
acc = np.mean(cross_val_score(clf, X_filtered, y, cv=5))
return acc
def show_relevance_redundancy(X, y, indices=None, title=""):
if not indices is None:
X = X.iloc[:, indices].copy()
y = y
fig = make_subplots(
rows=1,
cols=2,
column_widths=[0.68, 0.32],
column_titles=["relevance", "redundancy"],
)
trace_rel = px.bar(np.array([abs(np.corrcoef(x, y)[0, 1]) for x in X.values.T]))
trace_red = px.imshow(abs(np.corrcoef(X.values, rowvar=False)))
fig.add_trace(trace_rel.data[0], row=1, col=1)
fig.add_trace(trace_red.data[0], row=1, col=2)
fig.update_layout(width=1200, height=480, title=title)
fig.show()
show_relevance_redundancy(
X,
y,
None,
f"full dataset: acc={evaluate_model(clf, X, y, None):.4f}",
)
结果如下:
基线性能
模型的准确性为 0.9082。特征在相关性方面变化显著。有一些冗余的区域,特征之间似乎存在强相关性。
为了最小化目标函数,让我们创建一个带有约束的二次模型,并将其发送到量子退火器。代码将在下面解释。
k = 30
# Pearson correlation
correlation_matrix = abs(np.corrcoef(np.hstack((X, y[:, np.newaxis])), rowvar=False))
# fix the alpha parameter from the Milne paper
# to account for the numerous quadratic terms that are possible
beta = 0.5
alpha = 2 * beta * (k - 1) / (1 - beta + 2 * beta * (k - 1))
# generate weights for linear and quadratic terms, per Milne algorithm
Rxy = correlation_matrix[:, -1]
Q = correlation_matrix[:-1, :-1] * (1 - alpha)
np.fill_diagonal(Q, -Rxy * alpha)
# create binary quadratic model from the linear and quadratic weights
bqm = dimod.BinaryQuadraticModel(Q, "BINARY")
# create constrained quadratic model
cqm = dimod.ConstrainedQuadraticModel()
# the objective function of the CQM is the same as BQM
cqm.set_objective(bqm)
# constraint: limit the number of features to k
cqm.add_constraint_from_iterable(
((i, 1) for i in range(len(cqm.variables))), "==", rhs=k
)
# the sampler that will be used is the hybrid sampler in the DWave cloud
sampler = LeapHybridCQMSampler()
# solve the problem
sampleset = sampler.sample_cqm(cqm)
这里发生了很多事情。让我们解释每一步。
我们将特征数量限制为 k=30。这是模型的主要限制。
我们稍微偏离了论文中描述的目标函数。我们没有直接定义 alpha,而是使用了一个等效参数 beta,它具有相同的作用。然后我们以一种保持交互项贡献受控的方式定义 alpha——如果特征数量极大,这将确保交互项不会压倒目标函数中的线性项。
我们以一种可以直接构建 BQM 的方式塑造相关矩阵。我们将 BQM 限制在不超过 k=30 个特征的条件下,因此我们获得一个受限的二次模型 (CQM)。我们将 CQM 发送到量子退火器并收集结果。
应该注意的是,量子部分代码的运行时间约为 10 秒。对于许多问题,这是一个典型的基准时长,即使变量数量很大。运行时间往往随着问题复杂性的增加而增加,但这种增加可能不会像你从经典求解器中预期的那样急剧。
我们还没有完成。D-Wave 硬件通常会对解决方案空间进行采样,并返回大量可能的可行解决方案。这是因为硬件的工作方式:一个退火事件运行迅速且容易,因此重复进行退火是有意义的,这里自动发生了。如果生成足够多的样本,其中一些将是最佳解决方案。
所以我们需要对解决方案进行排序并挑选最佳方案。“最佳”意味着——它在目标函数中具有最低的值,D-Wave 称之为“能量”,因为它确实是量子处理器的物理能量。
# postprocess results
feasible = sampleset.filter(lambda s: s.is_feasible)
if feasible:
best = feasible.first
else:
assert len(cqm.constraints) == 1
best = sorted(
sampleset.data(),
key=lambda x: (list(cqm.violations(x.sample).values())[0], x.energy),
)[0]
assert list(best.sample.keys()) == sorted(best.sample.keys())
is_selected = np.array([bool(val) for val in best.sample.values()])
features = np.array([i for i, val in enumerate(is_selected) if val])
best_energy = best.energy
features
是一个包含量子退火器选定特征索引的数组。这是特征选择过程的结果。显然,它的长度将是 k=30。
让我们在特征选择后测量模型的准确性:
show_relevance_redundancy(
X,
y,
features,
f"explicit optimization: acc={evaluate_model(clf, X, y, features):.4f}",
)
显式优化
特征选择后的模型准确度为 0.9381。我们获得了确切数量的特征。大多数特征具有很高的相关性。特征之间的相关值通常较低。模型表现更好,可能更容易解释,并且选择正确特征所花费的时间并不长。
但是这个过程很长且繁琐,如果你只想选择一些特征的话。幸运的是,现在有一种更简单的方法。
轻松的特征选择方法
如果你安装了 D-Wave scikit-learn 插件,你只需执行以下操作:
X_new = SelectFromQuadraticModel(num_features=30, alpha=0.5).fit_transform(X, y)
这就是全部了。后台代码创建了二次模型,对其进行约束,发送到量子求解器,获取结果,再进行解析,最后选择最佳结果以 NumPy 数组格式返回,这是 scikit-learn 所期望的格式。运行时间大致相同。但让我们看看结果是什么样的:
X_new_df = pd.DataFrame(data=X_new, columns=list(range(X_new.shape[1])))
show_relevance_redundancy(
X_new_df,
y,
None,
f"plugin optimization: acc={evaluate_model(clf, X_new_df, y, None):.4f}",
)
插件优化
模型的准确度为 0.9369,本质上与我们通过显式优化获得的性能相同(在各种随机组件的典型波动范围内)。
选择的特征集略有不同。这可能是由于我的手动过程和库中的自动化实现之间的微小差异造成的。
在任何一种情况下,我们都将模型的性能从良好提升到了卓越,通过一种不进行启发式猜测的算法,而是考虑了所有特征的总和。
待探索的进一步主题
alpha 参数会使特征选择算法偏向于减少冗余但可能质量较低(alpha=0),或者偏向于提高质量但代价是增加冗余(alpha=1)。最佳值取决于你要解决的问题。
SelectFromQuadraticModel()
方法有一个名为 time_limit
的参数,正如名称所示:它控制求解器上的最大时间。你为所用的时间付费,因此这里的高值可能会很昂贵。另一方面,如果量子退火似乎未能收敛到足够好的解决方案,这里较高的值可能值得探索。这里展示的问题对量子处理器而言相当简单,所以我们仅使用了默认值。
请记住,你的计算机与云端量子求解器之间的数据交换可能并非简单,因此在时间敏感的应用程序(如金融)中,连接的延迟可能非常重要。
参考资料
D-Wave 网络研讨会介绍了 scikit-learn 插件。www.youtube.com/watch?v=VHEpe00AXPI
Milne, A., Rounds, M., & Goddard, P. (2018). 使用量子退火器进行信用评分和分类的最佳特征选择。1QBit.com。1QBit.com 上的白皮书
Boutell, M., Luo, J., Shen, X., Brown, C. (2004). 学习多标签场景分类。《模式识别》期刊,ScienceDirect。文章链接
入门指南:D-Wave 求解器
D-Wave scikit-learn 插件:D-Wave scikit-learn 插件
D-Wave 示例:CQM 特征选择:D-Wave 示例:CQM 特征选择
OpenML 场景识别数据集:OpenML 场景识别数据集
为本文创作的代码和工件:为本文创作的代码和工件
感谢罗斯-霍尔曼科技学院的马修·布特尔博士授权我在创作共用许可证下访问场景特征数据集。
所有图像均由作者创作。
使用 React 构建互动界面以展示令人兴奋的数据集
原文:
towardsdatascience.com/using-react-to-build-interactive-interfaces-to-exciting-dataset-c01691a5fc38
数据教程
使用网页开发创建更动态的数据可视化体验
·发表在 Towards Data Science ·5 分钟阅读·2023 年 9 月 19 日
–
界面演示视频
除了我作为一家小型机器学习公司的首席执行官的全职工作之外,我的爱好是创建美丽的数据可视化。
我通常使用 Matplotlib,但这次我想创造一个更互动的体验。
由于我喜欢网页开发和设计,我决定为世界银行的 人口估计和预测 数据集创建一个 React 应用程序。
这是一个令人着迷的数据集,你可以查看 1960 年至 2022 年所有国家和地区的人口金字塔,包括对 2050 年的预测。它的许可证是 Creative Commons Attribution 4.0。
这是一个非常适合互动界面的数据集,人们可以快速更改年份和地区。
在这个故事中,我将分享我的工作见解和所学到的经验。
如果你想测试这个解决方案,你可以在这里找到:datawonder.io/population-pyramids
让我们开始吧。
第一部分:准备数据
我想创建一个简单而快速的后台,将数据提供给前端,而不进行任何耗时的预处理。
相反,我的想法是提前处理所有数据,并在应用程序启动时将其全部加载到内存中。
世界银行的数据总是有一组指标,而我需要的指标具有以下格式:
人口年龄 ,
有 17 个年龄组,从 0–4 岁到 80 岁及以上。每个指标都有一个单独的列表示每一年,如下方的 pandas 数据框所示。
因为我确切知道我需要数据的哪些部分,不想在后台进行任何过滤或其他操作。
相反,我决定将 DataFrame 转换为以下格式的 JSON 文件。
data = {
<Area1> = {
"1960": {
"total": X + Y ... + A + B + ...
"male": {
"0-4": X,
"4-9": Y,
...
}
"female": {
"0-4": A,
"4-9": B,
...
}
}
...
}
...
}
这就是准备数据的全部内容。
第二部分:创建后端
后端是这个 web 应用程序中最简单的部分,只提供用户界面的 JSON 数据。
当应用程序启动时,我将预处理的 JSON 加载到内存中作为 Python 字典。
import json
from typing import List
from fastapi import FastAPI, Query
app = FastAPI()
population_data_file = open("./data/data.json")
population_data = json.load(population_data_file)
population_data_file.close()
我决定使用 FastAPI,因为它是一个可行的选项,并且使创建端点变得容易。
当用户发送查询时,我可以立即访问数据。
@app.get("/population")
def get_data(areas: List[str] = Query(None)):
populations = []
for area in areas:
if area not in population_data.keys():
return "Unknown area", 400
populations.append({"area": area, "data": population_data[area]})
return {"populations": populations}
我还有一个端点来列出所有可用的国家和地区。
@app.get("/areas")
def get_countries():
return list(population_data.keys())
这就是我后端所有的代码。
第三部分:创建用户界面
构建前端是这个项目中最耗时的部分。我不能详细讲解每一行代码。
我的目标是创建一个简约的应用程序,它既响应迅速又易于使用,因为许多数据接口过于复杂。
我使用了两个主要的 JavaScript 库来让我的生活更轻松:Recharts 和 Chakra。
我不会详细讲解 Recharts,因为它们有很棒的文档来描述你可以做的所有事情。
实质上,它们提供了一个叫做 ResponsiveContainer
的 React 组件,作为你图表的基础。
然后你可以定义像BarChart
、YAxis
、XAxis
、Tooltip
和Legend
这样的东西来创建完美的图表。
结构如下;所有组件都有几个属性来确保样式正确。
<ResponsiveContainer>
<BarChart>
<CartesianGrid />
<YAxis />
<XAxis />
<Tooltip />
<Legend />
<Bar />
<Bar />
</BarChart>
</ResponsiveContainer>
例如,我将我的 BarChart
组件设置为垂直布局,并将条形图堆叠在“sign”上,以使它们朝相反方向延展。
<BarChart
stackOffset="sign"
layout="vertical"
data={
prepare_data(
population.data[year.toString()],
percentage
)
}
>
<Bar dataKey="Male" stackId="a" fill="#57BED6" z={100} />
<Bar dataKey="Female" stackId="a" fill="#EC6237" z={10} />
</BarChart>
prepare_data()
函数将一个区域和年份的人口数据转换为对象列表。
const prepare_data = (population, percentage) => {
return AGES.map((age) => ({
name: age,
Male:
(population["male"][age] * -1) /
(percentage ? 0.01 * population["total"] : 1),
Female:
population["female"][age] / (percentage ? 0.01 * population["total"] : 1),
}));
};
下面是每个图表的样子。
我为其他组件设置了额外的属性来样式化轴标签和网格等。
第四部分:部署应用程序
我决定使用 Digital Ocean 的应用平台作为我的部署方式,但这可能不是最便宜的选择。
他们的最小后端服务比我需要的要好得多,但我计划随着时间的推移添加额外的端点以支持其他数据集。
部署前端
如果你使用create-react-app
创建前端,你可以将 Digital Ocean 指向那个仓库,它会在你将代码推送到 master 时构建应用程序。
你无需创建 Docker 容器或其他类似的东西。这对于像我这样在 DevOps 上比较懒的人来说非常棒。
它对于这个应用程序也很完美,因为它不需要任何花哨的东西。
结论
静态数据可视化很棒,但有时你会想要一种更互动的体验,让用户自己探索模式。
在这种情况下,构建简单的 React web 应用程序来使用开源库创建图表和图形是出奇的简单。
在本教程中,我分享了一些关于如何创建一个用户界面的见解,该界面允许用户查看不同国家和地区的年龄分布。
我希望你能从中获得一些灵感,自己做一个类似的项目。
感谢你的阅读,下次见!😊
使用无服务器函数来管理和监控基于云的训练实验
一种简单的常规方法,可以帮你节省大量金钱
·
查看 发表在 Towards Data Science · 11 min 阅读 · 2023 年 12 月 17 日
–
图片由 Ziyou Zhang 在 Unsplash 提供
这篇博客文章是与我的同事Shay Margalit共同撰写的。它总结了他关于如何利用AWS Lambda函数来提高对Amazon SageMaker训练服务的使用控制和成本的研究。感兴趣吗?请继续阅读 😃.
我们很幸运(或者说非常不幸——这要看你问谁)能亲历这场被许多人预期会改变我们所知世界的人工智能革命。得益于硬件开发的进步和对大量数据的获取,这场革命很可能会影响我们日常生活的许多方面——虽然具体如何,没有人能确定。为了满足对人工智能日益增长的需求,基础的机器学习模型规模迅速扩大,训练这些模型所需的资源也在增加。总而言之,要在人工智能开发领域保持相关性,就需要对昂贵的重型设备进行大规模投资。
基于云的托管训练服务,例如 Amazon SageMaker、Google Vertex AI 和 Microsoft Azure ML,通过使开发者能够在他们原本负担不起的机器上进行训练,降低了人工智能开发的门槛。虽然这些服务减少了人工智能的前期成本,并使你仅为训练时间付费,但变量成本的潜在累积需要对如何使用训练服务以及它们如何贡献于整体训练费用进行仔细规划。然而,不可避免的是,事情并不总是按计划进行。借用一句古老的意第绪谚语“开发者计划,编程之神笑”。当涉及到高风险时,比如训练人工智能模型——一个错误的实验可能导致数百或数千美元的计算时间浪费,因此明智之举是建立多个防线。
防线第一步——鼓励健康的发展习惯
第一防线应关注机器学习算法工程师的开发实践。以下是一些你可能会考虑的指导原则:
-
鼓励对用于训练的硬件资源进行适当且成本优化的使用(例如,参见这里)。
-
及早识别并终止失败的实验。
-
通过定期分析和优化运行时性能来提高性价比(例如,参见这里)。
尽管制定和调整像上面这些 AI 开发原则可能会提高您的生产力并减少浪费,但它们并不能完全防护所有可能的失败。例如,专用的故障检测运行时进程可能无法解决培训实验停滞的情况(例如,由于培训应用程序进程中的死锁),而培训任务仍然保持活动状态,直到它被主动停止或超时。
第二道防线 — 部署跨项目的保护措施
在这篇文章中,我们提议建立第二道防线,监控项目(或组织)中的所有培训活动,验证其是否符合预定的规则,并在发现异常培训实验时采取适当的措施。一种方法是使用专用的无服务器函数,这些函数会在培训任务的不同阶段触发,并编程以评估任务的状态,并根据需要停止或重新启动(可能会更改任务设置)任务。在接下来的章节中,我们将演示如何使用 AWS Lambda 作为对抗异常 Amazon SageMaker 培训实验的第二道防线的几个示例。
免责声明
尽管我们选择了 Amazon SageMaker 和 AWS Lambda 进行演示,但本文内容同样适用于其他服务,并且可以为它们实现类似的功能。请不要将我们选择这些服务解读为对其使用的推荐。云端培训有多种选项,每种选项都有其自身的优缺点。最适合您的选择将极大地取决于您项目的具体细节。
虽然我们将分享一些无服务器代码的 Python 示例,但我们不会详细讲解如何创建和部署它们作为 AWS Lambda 函数。与 AWS Lambda 交互的方式有很多种。我们建议读者参阅官方 AWS 文档以了解更多信息。
以下示例是为演示目的而创建的。它们可能需要修改以适应您项目的具体需求。在调整我们提议的解决方案之前,请确保充分理解所有代码细节及相关服务费用。重要的是,我们将分享的代码没有经过严格测试。任何涉及创建和调用多个 Lambda 函数以及Amazon CloudWatch 警报(如本文所述)的解决方案都需要进行适当的验证,以防止冗余/孤立工件的积累。
我们强烈建议您将本文的详细信息与最新的 AWS Lambda 文档和支持库的最新版本进行核对。
强制执行开发者合规性
虽然云治理对于成功和高效地使用云服务通常是至关重要的,但其实施有时可能具有挑战性。例如:Amazon SageMaker 包含一个用于向训练作业添加标签的 API。这些标签可以用来包含与 SageMaker 作业相关的元数据,比如训练项目的名称、开发阶段、当前试验的目标、运行作业的开发组或用户的名称等等。这些元数据可以用来收集统计信息,比如每个项目或组的开发成本。在下面的代码块中,我们展示了如何将多个标签应用于 SageMaker 训练作业:
from sagemaker.pytorch import PyTorch
tags = [{'Key': 'username', 'Value': 'johndoe'},
{'Key': 'model_name', 'Value': 'mnist'},
{'Key': 'training_phase', 'Value': 'finetune'},
{'Key': 'description', 'Value': 'fine tune final linear layer'}]
# define the training job with tags
estimator = PyTorch(
entry_point='train.py',
framework_version='2.1.0',
role='<arn role>',
py_version='py310',
job_name='demo',
instance_type='ml.g5.xlarge',
instance_count=1,
tags=tags
)
# deploy the job to the cloud
estimator.fit()
自然,这些标签只有在我们能够强制执行其应用时才有用。这就是AWS Lambda的作用。使用Amazon EventBridge,我们可以监控 SageMaker 训练作业状态的变化,并注册一个在每次变化时被触发的函数。在下面的代码块中,我们提出了一个 Python 例程,每次作业启动时都会验证特定的 SageMaker 标签。如果缺少某个标签,作业将被自动终止。事件的结构可以在这里找到。请注意使用(更详细的)SecondaryStatus字段来轮询训练作业的状态(而不是TrainingJobStatus)。
import boto3
def stop_training_job(training_job_name):
sm_client = boto3.client("sagemaker")
response = sm_client.stop_training_job(TrainingJobName=training_job_name)
assert response['ResponseMetadata']['HTTPStatusCode'] == 200
# TODO - optionally send an email notification
def enforce_required_tags(training_job_name, event):
event_tags = event['detail']['Tags']
if 'model_name' not in event_tags:
stop_training_job(training_job_name)
# define lambda handler
def sagemaker_event_handler(event, _):
job_name = event['detail']['TrainingJobName']
job_secondary_status = event['detail']['SecondaryStatus']
if job_secondary_status == 'Starting':
enforce_required_tags(job_name, event)
AWS 提供了多种创建 Lambda 函数的方法。有关详细信息,请参见AWS Lambda文档。创建后,请确保将该函数设置为EventBridge 规则的目标。
同样的功能可以用来强制执行额外的开发规则,这些规则旨在控制成本,比如:可以使用的实例类型、每个作业的最大实例数量、作业的最大运行时间等等。
停止停滞的实验
想象一下以下场景:你计划了一个大型的云端训练作业,将在八台每小时 $30 的 ML 计算实例上运行三天。为了完成这个任务,你已确保了 $17,280 的预算(8 个实例 x 每小时 $30 x 24 小时 x 3 天)。你在出发去度假三天的周末前启动了训练作业。当你从假期回来时,发现作业开始后一小时内,训练过程停滞,导致昂贵的设备实际上在三天内完全闲置。你不仅浪费了 $17,280(祝你好运向老板解释这一点),而且你的开发进度现在也被推迟了三天!!
保护自己免受这种情况的一种方法是监控底层训练作业资源的使用情况。例如,如果你的训练实例的 GPU 使用率在较长时间内低于某个阈值,这很可能是一个信号,表明出现了问题,需要立即停止训练作业。
我们将通过定义一个 Amazon CloudWatch 警报 来实现这一点,该警报监控每个 SageMaker 作业的一个训练实例的 GPU 使用情况,并在警报触发时调用一个 AWS Lambda 函数来终止该作业。设置这一功能需要三个组件:一个 Amazon CloudWatch 警报(每个训练作业一个),一个 AWS Lambda 函数,以及一个 Amazon Simple Notification Service (SNS) 主题,用于将 Lambda 函数与 CloudWatch 警报连接起来。
首先,我们创建一个 SNS 主题。这可以通过 Amazon SNS 控制台完成,也可以通过 Python 实现,如下所示:
import boto3
sns_client = boto3.client('sns')
# Create a SNS notification topic.
topic = sns_client.create_topic(Name="SageMakerTrainingJobIdleTopic")
topic_arn = topic.arn
print(f"Created SNS topic arn: {topic_arn}")
接下来,我们扩展了上面定义的 sagemaker_event_handler 函数,每次启动训练作业时创建一个独特的警报。我们编程设置警报以测量五分钟周期内的平均 GPU 使用率,并在连续三次测量低于 1% 时提醒我们的 SNS 主题。训练作业完成时,警报将被删除。
def create_training_alarm(job_name):
topic_arn = '<sns topic arn>'
SAMPLE_PERIOD_SECONDS = 60 * 5 # 5 minutes
SAMPLE_POINTS_LIMIT = 3
GPU_UTIL_THRESHOLD_PERCENTAGE = 1
cloudwatch_client = boto3.client('cloudwatch')
# A new sample is generated each SAMPLE_PERIOD_SECONDS seconds.
# The alarm will set off it there will be more than SAMPLE_POINTS_LIMIT
# below the limit.
response = cloudwatch_client.put_metric_alarm(
AlarmName=job_name + 'GPUUtil',
AlarmActions=topic_arn,
MetricName='GPUUtilization',
Namespace='/aws/sagemaker/TrainingJobs',
Statistic='Average',
Dimensions=[{
"Name": "Host",
"Value": job_name+"/algo-1"
}],
Period=SAMPLE_PERIOD_SECONDS,
EvaluationPeriods=SAMPLE_POINTS_LIMIT,
DatapointsToAlarm=SAMPLE_POINTS_LIMIT,
Threshold=GPU_UTIL_THRESHOLD_PERCENTAGE,
ComparisonOperator='LessThanOrEqualToThreshold',
TreatMissingData='notBreaching'
)
assert response['ResponseMetadata']['HTTPStatusCode'] == 200
def delete_training_alarm(job_name):
cloudwatch_client = boto3.client('cloudwatch')
response = cloudwatch_client.delete_alarms(
AlarmNames=[job_name+'GPUUtil'])
def sagemaker_event_handler(event, _):
job_name = event['detail']['TrainingJobName']
job_secondary_status = event['detail']['SecondaryStatus']
if job_secondary_status == 'Starting':
enforce_required_tags(job_name, event)
elif job_secondary_status == 'Training':
create_training_alarm(job_name)
elif job_secondary_status in ['Completed', 'Failed', 'Stopped']:
delete_training_alarm(job_name)
最后,我们定义一个第二个 Python AWS Lambda 函数,该函数 解析从 SNS 主题接收到的消息,并终止与警报关联的训练作业。
import boto3, json
def lambda_sns_handler(event, context):
data = json.loads(event['Records'][0]['Sns']['Message'])
alarm_name = data['AlarmName']
training_job_name = alarm_name.replace('GPUUtil', '')
stop_training_job(training_job_name)
AWS 提供了多种机制来订阅 Lambda 函数到 SNS 主题,包括 AWS 控制台、AWS CLI 和 AWS 无服务器应用程序模型 (AWS SAM)。
我们描述的解决方案总结在下图中:
AWS 架构图(作者提供)
请注意,相同的架构也可以用于强制执行你机器学习训练项目的最低 GPU 利用率。GPU 通常是训练基础设施中最昂贵的资源,你的目标应是最大化所有训练工作负载的利用率。通过规定最低利用率(例如 80%),你可以确保所有开发者适当优化他们的工作负载。
确保开发的连续性
在我们之前的示例中,我们演示了如何识别并停止一个停滞的实验。在我们描述的大型训练任务场景中,这有助于节省大量资金,但并没有解决开发的三天延迟。显然,如果停滞的源头在你的代码中,那么推迟恢复训练直到问题解决是有意义的。然而,我们经常遇到的训练中断并不是由于我们的代码造成的,而是由于服务环境中的偶发故障。在这种情况下,你的优先任务可能是确保训练的连续性,而不是等待有人手动恢复训练任务(使用最新的训练检查点)。在下面的代码块中,我们使用boto3的create_training_job API 来扩展我们的sagemaker_event_handler函数,以(简单地)恢复任何运行至少两小时后失败的训练任务。
import boto3, datetime
def clone_job(training_name, disable_spot=False):
# get description
client = boto3.client('sagemaker')
desc = client.describe_training_job(TrainingJobName=training_name)
# update the training name
new_training_name = training_name + 'clone'
use_spots = (not disable_spot) and desc["EnableManagedSpotTraining"]
if disable_spot:
desc["StoppingCondition"].pop("MaxWaitTimeInSeconds", None)
client.create_training_job(
TrainingJobName=new_training_name,
HyperParameters=desc["HyperParameters"],
AlgorithmSpecification=desc["AlgorithmSpecification"],
RoleArn=desc["RoleArn"],
OutputDataConfig=desc["OutputDataConfig"],
ResourceConfig=desc["ResourceConfig"],
StoppingCondition=desc["StoppingCondition"],
EnableNetworkIsolation=desc["EnableNetworkIsolation"],
EnableInterContainerTrafficEncryption=desc[
"EnableInterContainerTrafficEncryption"
],
EnableManagedSpotTraining=use_spots,
Tags=client.list_tags(ResourceArn=desc['TrainingJobArn'])
)
def sagemaker_event_handler(event, _):
TRAIN_TIME_THRESHOLD = 2 * 60 * 60: # 2 hours
job_name = event['detail']['TrainingJobName']
job_secondary_status = event['detail']['SecondaryStatus']
if job_secondary_status == 'Starting':
enforce_required_tags(job_name, event)
elif job_secondary_status == 'Training':
create_training_alarm(job_name)
elif job_secondary_status in ['Completed', 'Failed', 'Stopped']:
delete_training_alarm(job_name)
if job_secondary_status == 'Failed':
start_time = datetime.datetime.utcfromtimestamp(
event['detail']['CreationTime']/1000)
end_time = datetime.datetime.utcfromtimestamp(
event['detail']['TrainingEndTime']/1000)
training_time_seconds = (end_time - start_time).seconds
if training_time_seconds >= TRAIN_TIME_THRESHOLD:
clone_job(job_name)
上述函数会自动恢复任何在两小时后失败的任务。一个更实用的解决方案可能会尝试诊断错误类型,以确定恢复任务是否合适。一种方法是解析失败描述消息和/或与失败任务相关的 CloudWatch 日志。
高级 Spot 实例利用
Amazon SageMaker 的一个吸引人的特点是其对托管抢占式训练的支持。Amazon EC2 Spot 实例允许你利用未使用的 EC2 容量,并以折扣价格购买。这些实例的一个问题是它们在使用过程中可能会被收回(“中断”)。因此,Spot 实例应仅用于容错工作负载。SageMaker 通过代表你识别 Spot 中断,并在新的 Spot 实例可用时自动重新启动作业,使利用 Spot 实例变得简单。尽管托管抢占式 实例可以用来降低训练成本,但有时这一策略可能会适得其反。例如,当 Spot 容量不足时,你的训练作业可能会在开始之前超时。或者,作业可能会经历频繁的中断,这使得它无法取得任何有意义的进展。这两种情况都可能干扰开发并降低生产力。这些类型的情况可以使用 AWS Lambda 进行监控和处理。在下面的代码块中,我们扩展了我们的sagemaker_event_handler 函数,以识别一个被中断超过三次的训练作业,并用一个禁用了托管抢占式训练的克隆作业替换它。
def sagemaker_event_handler(event, _):
TRAIN_TIME_THRESHOLD = 2 * 60 * 60: # 2 hours
MIN_ITERRUPTS = 3
job_name = event['detail']['TrainingJobName']
job_secondary_status = event['detail']['SecondaryStatus']
if job_secondary_status == 'Starting':
enforce_required_tags(job_name, event)
elif job_secondary_status == 'Training':
create_training_alarm(job_name)
elif job_secondary_status in ['Completed', 'Failed', 'Stopped']:
delete_training_alarm(job_name)
if job_secondary_status == 'Failed':
start_time = datetime.datetime.utcfromtimestamp(
event['detail']['CreationTime']/1000)
end_time = datetime.datetime.utcfromtimestamp(
event['detail']['TrainingEndTime']/1000)
training_time_seconds = (end_time - start_time).seconds
if training_time_seconds >= TRAIN_TIME_THRESHOLD:
clone_job(job_name)
if job_secondary_status == 'Interrupted':
transitions = event['detail']["SecondaryStatusTransitions"]
interrupts = [e for e in transitions if e["Status"] == "Interrupted"]
num_interrupts = len(interrupts)
if num_interrupts > MIN_ITERRUPTS:
stop_training_job(job_name)
clone_job(job_name, disable_spot=True)
上述实现仅根据训练作业的中断次数确定了 Spot 使用策略。一个更复杂的解决方案可能还会考虑其他作业(使用相同实例类型)、中断发生的时间段、活跃训练时间的数量和/或由于 Spot 实例容量不足而超时的最近作业数量。
总结
有效的 AI 模型开发需要定义一个创意且详细的训练基础设施架构,以最小化成本并最大化生产力。在这篇文章中,我们展示了如何使用无服务器的 AWS Lambda 函数来增强 Amazon SageMaker 的托管训练服务,以解决训练过程中可能出现的一些常见问题。当然,具体应用这些技术的方式将很大程度上取决于你项目的具体情况。
如果有任何问题、意见或修正,请随时联系。务必查看我们关于 DL 训练优化的其他文章。
使用 SHAP 调试 PyTorch 图像回归模型
原文:
towardsdatascience.com/using-shap-to-debug-a-pytorch-image-regression-model-4b562ddef30d
使用 DeepShap 来理解和改进驱动自动驾驶汽车的模型
·发表于 Towards Data Science ·阅读时间 11 分钟·2023 年 1 月 10 日
–
(来源:作者)
自动驾驶汽车让我感到恐惧。巨大的金属块在空中飞驰,如果出现问题没有人可以阻止它们。为了降低风险,仅仅评估这些“怪物”背后的模型是不够的。我们还需要了解它们是如何进行预测的。这是为了避免任何可能导致意外的边缘情况。
好吧,我们的应用程序并没有那么重要。我们将调试用于驱动迷你自动驾驶汽车的模型(最坏的情况可能是一个瘀伤的脚踝)。尽管如此,IML 方法仍然有用。我们将看看它们如何甚至能提高模型的性能。
具体来说,我们将:
-
使用 PyTorch 和图像数据以及连续目标变量来微调 ResNet-18
-
使用 MSE 和散点图来评估模型
-
使用 DeepSHAP 解释模型
-
通过更好的数据收集来纠正模型
-
讨论图像增强如何进一步改善模型
在这个过程中,我们将讨论一些关键的 Python 代码。你还可以在GitHub上找到完整的项目。
如果你是 SHAP 的新手,请查看下面的视频。如果你想要更多内容,可以查看我的SHAP 课程。 如果你注册我的新闻通讯,你可以免费访问 😃
Python 包
# Imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import glob
import random
from PIL import Image
import cv2
import torch
import torchvision
from torchvision import transforms
from torch.utils.data import DataLoader
import shap
from sklearn.metrics import mean_squared_error
数据集
我们从在一个房间里收集数据开始这个项目(这将对我们造成困扰)。如前所述,我们使用图像来驱动自动化汽车。你可以在Kaggle上找到这些图像的例子。这些图像的尺寸均为 224 x 224 像素。
我们用下面的代码显示了其中一个。注意图像名称(第 2 行)。前两个数字是 224 x 224 帧中的 x 和 y 坐标。在图 1中,你可以看到我们使用绿色圆圈显示了这些坐标(第 8 行)。
#Load example image
name = "32_50_c78164b4-40d2-11ed-a47b-a46bb6070c92.jpg"
x = int(name.split("_")[0])
y = int(name.split("_")[1])
img = Image.open("../data/room_1/" + name)
img = np.array(img)
cv2.circle(img, (x, y), 8, (0, 255, 0), 3)
plt.imshow(img)
图 1:赛道输入图像示例(来源:作者)
这些坐标是目标变量。模型使用图像作为输入来预测它们。然后,这一预测用于引导汽车。在这种情况下,你可以看到汽车正在接近左转。理想的方向是朝向绿色圆圈所给的坐标。
训练 PyTorch 模型
我想专注于 SHAP,因此我们不会深入讨论建模代码。如果你有任何问题,随时在评论中提问。
我们首先创建ImageDataset类。这个类用于加载我们的图像数据和目标变量。它通过paths来完成这一任务。值得指出的是目标变量的缩放方式——x和y的范围都在**-1和1**之间。
class ImageDataset(torch.utils.data.Dataset):
def __init__(self, paths, transform):
self.transform = transform
self.paths = paths
def __getitem__(self, idx):
"""Get image and target (x, y) coordinates"""
# Read image
path = self.paths[idx]
image = cv2.imread(path, cv2.IMREAD_COLOR)
image = Image.fromarray(image)
# Transform image
image = self.transform(image)
# Get target
target = self.get_target(path)
target = torch.Tensor(target)
return image, target
def get_target(self,path):
"""Get the target (x, y) coordinates from path"""
name = os.path.basename(path)
items = name.split('_')
x = items[0]
y = items[1]
# Scale between -1 and 1
x = 2.0 * (int(x)/ 224 - 0.5) # -1 left, +1 right
y = 2.0 * (int(y) / 244 -0.5)# -1 top, +1 bottom
return [x, y]
def __len__(self):
return len(self.paths)
实际上,当模型部署时,仅使用 x 预测来引导汽车。由于缩放,x 预测的符号将决定汽车的方向。当x < 0时,汽车应该左转。同样,当x > 0时,汽车应该右转。x 值越大,转弯越急。
我们使用 ImageDataset 类创建训练和验证数据加载器。这是通过对所有来自房间 1 的图像路径进行随机80/20拆分来完成的。最终,我们在训练集和验证集中分别有1,217和305张图像。
TRANSFORMS = transforms.Compose([
transforms.ColorJitter(0.2, 0.2, 0.2, 0.2),
transforms.Resize((224, 224)),
transforms.ToTensor(),
transforms.Normalize([0.485, 0.456, 0.406], [0.229, 0.224, 0.225])
])
paths = glob.glob('../data/room_1/*')
# Shuffle the paths
random.shuffle(paths)
# Create a datasets for training and validation
split = int(0.8 * len(paths))
train_data = ImageDataset(paths[:split], TRANSFORMS)
valid_data = ImageDataset(paths[split:], TRANSFORMS)
# Prepare data for Pytorch model
train_loader = DataLoader(train_data, batch_size=32, shuffle=True)
valid_loader = DataLoader(valid_data, batch_size=valid_data.__len__())
注意valid_loader的batch_size。我们使用验证数据集的长度(即 305)。这使我们可以在一次迭代中加载所有验证数据。如果你处理的是更大的数据集,可能需要使用较小的批量大小。
我们加载一个预训练的 ResNet18 模型(第 5 行)。通过设置model.fc,我们更新了最后一层(第 6 行)。这是一个从 512 个节点到我们 2 个目标变量节点的全连接层。我们将使用 Adam 优化器来微调这个模型(第 9 行)。
output_dim = 2 # x, y
device = torch.device('mps') # or 'cuda' if you have a GPU
# RESNET 18
model = torchvision.models.resnet18(pretrained=True)
model.fc = torch.nn.Linear(512, output_dim)
model = model.to(device)
optimizer = torch.optim.Adam(model.parameters())
我已经使用 GPU 训练了模型(第 2 行)。你仍然可以在 CPU 上运行代码。微调不像从头开始训练那样计算密集!
最后,我们有了模型训练代码。我们使用 MSE 作为损失函数训练 10 个周期。我们的最终模型是验证集上 MSE 最低的模型。
name = "direction_model_1" # Change this to save a new model
# Train the model
min_loss = np.inf
for epoch in range(10):
model = model.train()
for images, target in iter(train_loader):
images = images.to(device)
target = target.to(device)
# Zero gradients of parameters
optimizer.zero_grad()
# Execute model to get outputs
output = model(images)
# Calculate loss
loss = torch.nn.functional.mse_loss(output, target)
# Run backpropogation to accumulate gradients
loss.backward()
# Update model parameters
optimizer.step()
# Calculate validation loss
model = model.eval()
images, target = next(iter(valid_loader))
images = images.to(device)
target = target.to(device)
output = model(images)
valid_loss = torch.nn.functional.mse_loss(output, target)
print("Epoch: {}, Validation Loss: {}".format(epoch, valid_loss.item()))
if valid_loss < min_loss:
print("Saving model")
torch.save(model, '../models/{}.pth'.format(name))
min_loss = valid_loss
评估指标
此时,我们想了解我们的模型表现如何。我们查看 MSE 和实际与预测 x 值的散点图。暂时忽略 y,因为它不会影响汽车的方向。
训练和验证集
图 2提供了训练和验证集的这些指标。对角线红线表示完美预测。对于x < 0和x > 0,此线周围有类似的变异。换句话说,模型能够以类似的准确性预测左转和右转。在训练和验证集上的类似表现也表明模型没有过拟合。
图 2:模型在训练和验证集上的评估(来源:作者)
要创建上述图,我们使用model_evaluation函数。注意,数据加载器应创建为在第一次迭代中加载所有数据。
def model_evaluation(loaders,labels,save_path = None):
"""Evaluate direction models with mse and scatter plots
loaders: list of data loaders
labels: list of labels for plot title"""
n = len(loaders)
fig, axs = plt.subplots(1, n, figsize=(7*n, 6))
# Evalution metrics
for i, loader in enumerate(loaders):
# Load all data
images, target = next(iter(loader))
images = images.to(device)
target = target.to(device)
output=model(images)
# Get x predictions
x_pred=output.detach().cpu().numpy()[:,0]
x_target=target.cpu().numpy()[:,0]
# Calculate MSE
mse = mean_squared_error(x_target, x_pred)
# Plot predcitons
axs[i].scatter(x_target,x_pred)
axs[i].plot([-1, 1],
[-1, 1],
color='r',
linestyle='-',
linewidth=2)
axs[i].set_ylabel('Predicted x', size =15)
axs[i].set_xlabel('Actual x', size =15)
axs[i].set_title("{0} MSE: {1:.4f}".format(labels[i], mse),size = 18)
if save_path != None:
fig.savefig(save_path)
使用下面的函数可以看出我们的意思。我们创建了一个新的train_loader,将批量大小设置为训练数据集的长度。加载保存的模型也很重要(第 2 行)。否则,你将使用上一个纪元训练的模型。
# Load saved model
model = torch.load('../models/direction_model_1.pth')
model.eval()
model.to(device)
# Create new loader for all data
train_loader = DataLoader(train_data, batch_size=train_data.__len__())
# Evaluate model on training and validation set
loaders = [train_loader,valid_loader]
labels = ["Train","Validation"]
# Evaluate on training and validation set
model_evaluation(loaders,labels)
移动到新位置
结果看起来不错!我们预计汽车会表现良好,它确实如此。直到我们将它移动到新位置:
图 3:模型在新位置表现不佳(来源:作者)
我们从新位置(房间 2 和房间 3)收集了一些数据。对这些图像进行评估时,你会发现我们的模型表现得不如预期。这很奇怪!汽车在完全相同的轨道上,为什么房间会有影响呢?
图 3:模型在房间 2 和房间 3 上的评估(来源:作者)
使用 SHAP 调试模型
我们寻求 SHAP 的答案。它可以用来理解哪些像素对给定的预测重要。我们首先加载保存的模型(第 2 行)。SHAP 尚未实现 GPU 支持,所以我们将设备设置为 CPU(第 5-6 行)。
# Load saved model
model = torch.load('../models/direction_model_1.pth')
# Use CPU
device = torch.device('cpu')
model = model.to(device)
计算 SHAP 值时,我们需要获取一些背景图像。SHAP 将在计算值时对这些图像进行集成。我们使用batch_size为 100 的图像。这应该能给我们合理的近似值。增加图像数量将提高近似精度,但也会增加计算时间。
#Load 100 images for background
shap_loader = DataLoader(train_data, batch_size=100, shuffle=True)
background, _ = next(iter(shap_loader))
background = background.to(device)
我们通过将模型和背景图像传入DeepExplainer函数来创建一个解释器对象。这个函数有效地为神经网络近似 SHAP 值。作为替代,你可以用GradientExplainer函数替换它。
#Create SHAP explainer
explainer = shap.DeepExplainer(model, background)
我们加载 2 张示例图像——一张右转和一张左转(第 2 行),并进行变换(第 6 行)。这很重要,因为图像应该与训练模型时使用的格式相同。然后,我们计算这些图像预测的 SHAP 值(第 10 行)。
# Load test images of right and left turn
paths = glob.glob('../data/room_1/*')
test_images = [Image.open(paths[0]), Image.open(paths[3])]
test_images = np.array(test_images)
test_input = [TRANSFORMS(img) for img in test_images]
test_input = torch.stack(test_input).to(device)
# Get SHAP values
shap_values = explainer.shap_values(test_input)
最后,我们可以使用image_plot函数来显示 SHAP 值。但我们首先需要重新构造这些值。SHAP 值的返回维度是:
( #targets, #images, #channels, #width, #height)
我们使用转置函数,所以我们有维度:
(#targets, #images, #width, #height, #channels)
请注意,我们也将原始图像传入了image_plot函数。由于变换,test_input图像可能会显得很奇怪。
# Reshape shap values and images for plotting
shap_numpy = list(np.array(shap_values).transpose(0,1,3,4,2))
test_numpy = np.array([np.array(img) for img in test_images])
shap.image_plot(shap_numpy, test_numpy,show=False)
你可以在图 4中看到结果。第一列给出原始图像。第二列和第三列分别是 x 和 y 预测的 SHAP 值。蓝色像素降低了预测值。相比之下,红色像素增加了预测值。换句话说,对于 x 预测,红色像素导致了更尖锐的右转。
图 4:左转和右转的示例 shap 值(来源:作者)
现在我们有了进展。重要的结果是模型正在使用背景像素。你可以在图 5 中看到这一点,我们对右转的 x 预测进行了放大。换句话说,背景对预测很重要。这解释了表现不佳的原因!当我们转到一个新房间时,背景发生了变化,我们的预测变得不可靠。
图 5:右转 x 预测的 shap 值(来源:作者)
模型对房间 1 的数据过拟合。每张图像中都有相同的对象和背景。因此,模型将这些与左转和右转关联起来。由于训练和验证图像中都有相同的背景,我们在评估中无法识别出这一点。
图 6:对训练数据的过拟合(来源:作者)
改进模型
我们希望我们的模型在所有条件下表现良好。为此,我们期望它只使用轨迹上的像素。那么,让我们讨论一些提高模型鲁棒性的方法。
收集新数据
最好的解决方案是简单地收集更多的数据。我们已经有了一些来自房间 2 和 3 的数据。按照相同的过程,我们使用来自所有 3 个房间的数据训练一个新模型。查看图 7,它现在在新房间的图像上表现更好。
图 7:在房间 2 和 3 上评估新模型(来源:作者)
希望通过在多个房间的数据上训练,我们可以打破转弯与背景之间的关联。现在,左转和右转上存在不同的对象,但轨迹保持不变。模型应该学会轨迹才是预测的重要因素。
我们可以通过查看新模型的 SHAP 值来确认这一点。这些值对应于我们在图 4中看到的相同转弯。现在,背景像素的权重较少。好吧,虽然不完美,但我们在进步。
图 8:在所有 3 个房间上训练的模型的 shap 值(来源:作者)
我们可以继续收集数据。我们收集的数据地点越多,我们的模型就会越强大。然而,数据收集可能是耗时的(而且无聊!)。相反,我们可以考虑数据增强。
数据增强
数据增强是指我们使用代码系统地或随机地改变图像。这使我们能够人为地引入噪声并增加数据集的大小。
例如,我们可以通过翻转图像在垂直轴上来将数据集的大小翻倍。我们之所以可以这样做,是因为我们的轨道是对称的。如图 9 所示,删除也可能是一个有用的方法。这涉及到包含那些对象或整个背景被去除的图像。
图 9:使用删除进行图像增强的示例(来源:作者)
在构建强大的模型时,您还应该考虑诸如光照条件和图像质量等因素。我们可以通过颜色抖动或添加噪声来模拟这些因素。如果您想了解所有这些方法,请查看下面的文章。
使用 Python 进行数据增强,包括翻转、调整亮度、颜色抖动和随机噪声
towardsdatascience.com
在上述文章中,我们还讨论了为什么很难判断这些增强是否使模型更强大。我们可以在许多环境中部署模型,但这很耗时。幸运的是,SHAP 可以作为一种替代方案。与数据收集一样,它可以帮助我们了解增强如何改变模型的预测方式。
希望您喜欢这篇文章!您可以通过成为我的推荐会员 😃 来支持我。
[## 通过我的推荐链接加入 Medium — Conor O’Sullivan
作为 Medium 会员,您的部分会员费会分配给您阅读的作者,同时您可以完全访问每一个故事……
conorosullyds.medium.com](https://conorosullyds.medium.com/membership?source=post_page-----4b562ddef30d--------------------------------)
| Twitter | YouTube | Newsletter — 免费注册以获得Python SHAP 课程
数据集
JatRacer 图像(CC0: 公共领域) www.kaggle.com/datasets/conorsully1/jatracer-images
参考资料
SHAP,PyTorch 深度解释器 MNIST 示例 shap.readthedocs.io/en/latest/example_notebooks/image_examples/image_classification/PyTorch%20Deep%20Explainer%20MNIST%20example.html
使用斜率图表简化你的数据可视化
原文:
towardsdatascience.com/using-slope-charts-to-simplify-your-data-visualization-be1f0eaf1f0f
数据可视化,数据讲述
通过使用斜率图表来简化你繁杂的图表:一个 Python Altair 教程
·发布于 Towards Data Science ·阅读时间 5 分钟·2023 年 12 月 8 日
–
作者提供的图片
我们可能会绘制图表以包含尽可能多的概念在我们的可视化中。因此,我们的图表可能会很难阅读且具有干扰性。因此,在绘制任何内容之前,请坐下来计划你想要传达的信息。然后,查看你的数据,决定什么是有效的必要内容。将其余部分排除在你的可视化之外。
在本教程中,我们将看看如何使用斜率图表来简化一个繁杂的趋势线。如果你是数据分析师,你可能会因为使用斜率图表看到显著的信息损失而惊慌。但我向你保证,在某些情况下,这确实是值得的。
让我们看看斜率图表可以应用的场景。
你可以在何时使用斜率图表
斜率图表是一种只显示第一个和最后一个点的折线图,如下面的图示。
作者提供的图片
当你只想了解数据趋势的斜率时,斜率图表特别有用。因此,斜率图表有助于简化趋势线。 例如,你可以使用斜率图表来查看产品销售在一段时间内是增加还是减少。假设你有许多趋势线需要在同一图表中表示,而你只对每条趋势线的第一个和最后一个值感兴趣。你可以通过使用斜率图表来简化图表。
让我们实现一个实际的例子,看看如何在 Python 数据可视化库 Altair 中实现斜率图表。
示例:游客到达情况
考虑旅游住宿设施的到达数据集,这是由 Eurostat 发布的开放数据。假设你想将葡萄牙的游客到达情况与其他五个国家进行比较:德国、法国、意大利、英国和西班牙。
从将数据集导入为 Pandas DataFrame 开始,并使用melt()
函数来解压数据集,即将数据行转换为列。
import altair as alt
import pandas as pd
from ydata_profiling import ProfileReport
df = pd.read_csv('tourist_arrivals_countries.csv', parse_dates=['Date'])
df = pd.melt(df, id_vars='Date', value_name='Tourist Arrivals', var_name='Country')
接下来,使用 Python Altair 绘制基本图表。要了解 Altair 的介绍,请阅读这篇介绍文章由Soner Yıldırım提供。
chart = alt.Chart(df).mark_line().encode(
x = 'Date:T', #A
y = 'Tourist Arrivals:Q', #B
color=alt.Color('Country:N') #C
)
chart.save('raw-chart.html')
下图显示了结果图表:
图片由作者提供
图表难以阅读。因此,我们通过以下方法改进它:
-
去除缺失的年份: 过滤掉 1994 年之前和 2018 年之后的所有年份
-
按年分组并计算总和: 不是显示所有月份的细节,而是仅显示年度值。
以下代码片段展示了如何进行操作:
# extract year from date
df.loc[:, 'Year'] = df['Date'].dt.year
# filter out years before 1994 and after 2018
df = df[(df['Year'] >= 1994) & (df['Year'] <= 2018)]
# group by year and country
df = df.groupby(['Year', 'Country'])['Tourist Arrivals'].sum().reset_index()
现在,再次绘制图表:
chart = alt.Chart(df).mark_line().encode(
x = 'Year:O',
y = 'Tourist Arrivals:Q',
color=alt.Color('Country:N')
)
chart.save(‘chart.html’)
下图显示了结果图表:
图片由作者提供
从图表中,我们注意到所有国家的趋势线都有所上升。因此,如果我们想评估趋势线在最后一年(2018 年)和第一年(1984 年)之间的变化,我们可以排除所有中间年份(1985 年至 2017 年),仅显示第一个和最后一个值。
此外,为了更清楚地展示差距,我们可以计算 2018 年与 1984 年之间的百分比增长,如下代码所示:
# select only 1994 and 2018
df = df[(df['Year'] == 1994) | (df['Year'] == 2018)]
# add a new column containing the difference for each country between the number of tourist arrivals in the year and 1994
for country in df['Country'].unique():
current = df[df['Country'] == country]['Tourist Arrivals']
base = df[(df['Country'] == country) & (df['Year'] == 1994)]['Tourist Arrivals'].values[0]
df.loc[df['Country'] == country, 'PI'] = (current - base)/ base*100
现在,让我们构建图表。我们还通过在strokeWidth
通道中使用条件来突出显示与葡萄牙相关的趋势线。
base = alt.Chart(df).encode(
x = alt.X('Year:O', title=''),
y = alt.Y('Tourist Arrivals:Q', title='Tourist Arrivals'),
color=alt.Color('Country:N',
scale=alt.Scale(scheme='set2'),
legend=None),
strokeWidth=alt.condition(alt.datum.Country == 'PT', alt.value(7), alt.value(0.5))
).properties(
title='Tourist Arrivals in Portugal (1994-2018)',
width=600,
)
chart = base.mark_line()
chart.save('chart.html')
下图显示了结果图表:
图片由作者提供
只需包括一个要点:与每个图表相关的标签。让我们将它们作为文本注释添加,如下代码片段所示:
annotation = base.mark_text(
dx=10,
align='left',
fontSize=12
).encode(
# format the text to show only 2 decimal places and add a percentage sign
text='label:N'
).transform_filter(
alt.datum.Year == 2018
).transform_calculate(
label='datum.Country + "(Increase:" + format(datum.PI, ".2f") + "%)"'
)
chart = (chart + annotation)
使用transform_calculate()
来正确格式化标签。下图显示了结果图表:
图片由作者提供
你可以继续优化图表,例如旋转 x 轴标签。不管怎样,斜率图已经完成,非常清晰 😃
总结
恭喜你!你刚刚学习了如何在 Python Altair 中实现斜率图!
使用斜率图仅显示趋势线的第一个和最后一个值。
你可以在我的书《数据讲故事与生成 AI:使用 Python 和 Altair》的GitHub 仓库中找到描述场景的完整代码。此外,这个例子是书中一个练习的解决方案。敬请关注,以实现更多示例 😃
如果你已经读到这里,今天的内容我觉得非常满意。谢谢,下次见 😃
附加资源
[书籍] 数据讲故事与生成 AI:使用 Python 和 Altair
你可能也会感兴趣…
这篇文章将展示如何使用 Python Altair 库构建地理地图。
在 Python Altair 中构建地理地图的 3 种方法
对于数据科学家来说,数据可视化是一项基本技能。它有助于快速理解数据中的模式和关联,否则可能会被忽视。地理地图是一个很好的方式来 继续阅读…
你可能不知道的 3 种时间序列可视化方法
在这篇文章中,我们将描述三种时间序列可视化的替代方法:
-
日历热力图
-
箱形图
-
循环图 继续阅读…
使用符号回归为 Elo 著名评分系统增加不确定性
并创建一个意外有用的评分算法
·发表于Towards Data Science ·阅读时间 12 分钟·2023 年 7 月 6 日
–
照片由JESHOOTS.COM提供,来自Unsplash
一个通用评分系统
Elo 评分系统在一些领域变得非常著名。最著名的例子可能是它自 1960 年代以来作为国际象棋评分的基础。此外,网站 538 成功地对其进行了修改,用于他们大部分知名的体育评分。较少公开的是,许多视频游戏开发者在他们的匹配系统中使用了 Elo 系统的变种。如果你正在阅读这篇文章,我会假设你对该系统有一定的了解。它为何在如此多的领域被广泛使用?我认为是因为它的计算扩展性、多功能性和简单性。然而,它也有一些缺点。本文将解决一个非常关键的问题,同时保持上述优点。
符号回归
尽管大型语言模型目前获得了所有的关注(双关含义),但也有其他令人兴奋的模型在被单独开发,具有非常不同的用途。符号回归通常适合发现封闭形式的解析规则,而不是处理像图像分类或音频翻译这样的深度学习任务。例如,如果你想重新发现牛顿的冷却定律,你可以构建一个资源密集型的密集神经网络。虽然这样做在数据足够的情况下效果良好,但它无法推广到未见过的情况。然而,符号回归将是正确的工具。它可以用有限的数据找到精确的公式,因此不仅可以推广,还能节省相当多的计算资源。Cranmer 等人撰写的我最喜欢的论文之一对此进行了更深入的探讨,甚至发展了一个以前未发现的暗物质过密度方程。
问题
经典的 Elo 评分将每个评分视为同样确定。这通常对大规模评分系统而言是一个不准确的假设。简单来说,新进入评分系统的玩家应该几乎总是以比那些已存在一段时间的玩家更大的方差进行建模。同样,那些评分系统长时间未曾见过的玩家可能已经有所改进或退步。假设你和你的四个密友经常在马里奥卡丁车中互相竞技。现在,我们可以使用简化的假设,你们只进行一对一的比赛。虽然还没有进行任何比赛,但假设你们有以下基本真实的 Elo 评分:
-
爱丽丝:1900
-
鲍勃:1700
-
切尔西:1500
-
德米特里:1300
-
伊夫琳:1100
在经典的 Elo 系统中,预期的胜率完全由评分差异决定:
其中 WP_A 是玩家 A 的胜率,R_B 和 R_A 是每个玩家的评分
因为这是一个模拟,我们不必担心当前的状态、良好的睡眠或血液酒精含量。相反,我们真实的 Elo 评分通过上述方程对应真实的胜率。这将意味着在模拟中每个玩家对抗另一名玩家时,以下是真实的胜率百分比:
你决定使用 Elo 系统来确定每个人的水平,因为目前你不知道真实评分。Elo 系统所需的唯一参数是 K 值。K 值是每场比赛中玩家之间“下注”的评分点数,最终根据赛前获胜概率奖励给胜者。较高的 K 值用于更快响应近期结果,而较低的 K 值用于更信任现有评分。在大多数使用情况下,你会看到 K 值在 10 到 40 之间。高于 40,评分变得非常嘈杂。低于 10,评分几乎不会变化。以下是 5,000 次模拟中,你平均需要多长时间才能确定每个人的水平,使用 K 值为 12、24 和 36,经过大约 100 场比赛:
如你所见,即使使用相当高的 K 值(36),我们的系统仍需 100 场比赛才能收敛。随着 K 值的增加,一致性变得更差。例如,当 K=12 时,我们对绿色玩家 Chelsea 的评分近乎真实评分 1500。然而,使用较高的 K 值时,我们对 Chelsea 评分的估计会在 1400 到 1600 之间漂移。
一个选择是从高 K 值开始,以提高收敛速度,然后随着评分的准确性提高而减少 K 值。这是一种在更多数据到来时减少不确定性的巧妙方法,但结果显示,它的表现不如更有原则的方法。
其他系统
过去曾有很多尝试解决这个问题,特别是 1995 年的Glicko 系统和 2005 年的TrueSkill 系统。Glicko 依赖于一个封闭形式的更新方程,而 TrueSkill 使用贝叶斯因子图来更新评分。TrueSkill 具有处理多人比赛场景或团队场景的巨大优势,并且通常收敛更快。TrueSkill 通常需要额外的计算,尤其是当每场比赛/比赛/锦标赛的玩家数量较多时。它也比 Glicko 对参数变化更敏感。使用哪种系统的决定当然取决于应用(但通常最终选择 TrueSkill)。让我们看看 Glicko、TrueSkill、Elo 以及带有衰减 K 的 Elo 如何尝试收敛到 Alice 的评分:
如你所见,Glicko 和 TrueSkill 在平均情况下收敛快得多,比 Elo 和带有 K 值衰减的 Elo 更快。从根本上说,更快的收敛是因为它们融入了不确定性度量。TrueSkill 在我的模拟中,实际超出真实评分,这显示了该系统对参数选择的敏感性。如果我调整参数,它将优于 Glicko,并且不会超出真实评分。
符号回归与我自己的系统
在我了解到其他人如何解决这个问题之前,我决定自己解决这个问题。我很高兴没有看到他们的解决方案,因为这可能会使我放弃创造自己的方法。在这个过程中,我发现了一个很好的使用案例,就是常常被低估的符号回归。
入门数学
与任何经典模型一样,我们必须从假设开始。就本文而言,可以安全地假设玩家和团队的技能在大多数情况下是正态分布的。根据我的经验,这对于几乎任何竞争都是接近真实的。国际象棋评级被发现略微更适合逻辑分布,但我现在想避免过于详细。
假设我们有一个正态分布的 10,000 名马里奥卡丁车玩家的样本。我们来放大其中的两个玩家。其中一个是“蓝色玩家”,另一个是“红色玩家”。
我会说,在这种情况下,我的表现比较符合红色玩家(从差到更差)。
在这种情况下,我从这个分布中随机抽取了两个玩家。在经典的 Elo 公式中,比赛完全可以通过这两个玩家之间的评分差异来描述。经典 Elo 没有试图考虑我们可能对蓝色玩家的数据比红色玩家少。为了澄清,我们还没有使用这些分布来预测单场比赛的结果。这些分布仅仅代表了整体的技能水平。即使它们完全不重叠,红色玩家击败蓝色玩家的概率仍然很小。
比如说,我们只见过蓝色玩家一次,他或她表现得相当不错。可能是因为在正确的时间拾取了正确的道具箱,这很难说。另一方面,我们对红色玩家的表现很确定,虽然具体有多差还有些不确定。在这种新的 Elo 形式中,我们承认蓝色玩家实际上可能比红色玩家差(因为分布有些有意义地重叠)。
让我们把红色分布叫做 R,蓝色分布叫做 B,差异分布叫做 Z:
我保证,数学方程不多,而且只有简单的方程
由于这些是独立分布,B 和 R 的协方差为零。我们的方程此时变得非常简单。换句话说,经典 Elo 只关心均值 Z,约为 1.85。在这种情况下,我们还跟踪 Z 的标准差,约为 0.84。这样,当我们记录一个观察值时,可以以贝叶斯方式更新分布。
注意上面的 x 轴具有相同的刻度。我只是通过 sigmoid 函数将技能差异转化为获胜概率。
上面,我可视化了当我们从技能差异分布中运行 Elo 10,000 次时的情况。经典 Elo 会给蓝色玩家一个 86.4%的胜率,用黄色星星表示。我们的新 Elo 模型将根据从分布中抽样的结果给出许多不同的概率。有趣的是,由于一个尾部场景使胜率急剧下降,而另一个尾部场景达到了上限,因此我们分布的平均胜率实际上比仅仅一个点估计要低得多——83.7%。我只是猜测,但这也许有助于解释为什么国际象棋联合会发现弱者往往比 Elo 预测的赢得更多。
还值得快速提及的是,一些游戏比其他游戏更具运气成分。例如,你可以在 Uno、卡坦岛、尤克或者 Five Crowns 中变得更出色,但这些游戏中涉及的运气因素足够多,以至于一个差的玩家仍然可能击败排名更高的玩家。另一方面,如果我与 Magnus Carlson 下棋,我 100 场 100 次都会输,因此涉及的运气成分较少。为了建模不可减少的运气因素,需要额外的复杂性(我现在暂时避免)。
符号回归
现在,我们如何利用胜负结果来更新这些分布的参数呢?我们需要将其作为一个二元结果变量,以保持 Elo 系统的通用性——我们希望能够将这个系统应用于几乎所有的游戏(一个好的练习是创建一个基于连续结果更新的系统)。然而,据我所知,没有简单的闭式解来更新带有二元结果的正态分布的参数,因此会出现一些不完美之处。一方面,我可以花费几年时间通过大学课程来磨练我的数学技能,以至于知道如何坐下来解决这个问题。另一方面,我可以作弊,把它放入符号回归模型中,并额外写一篇 Medium 文章。我选择了后者,尽管我非常尊重那些选择前者的人。
Beta 分布非常适合建模二元结果的概率,我们已经通过对正态分布均值的 sigmoid 变换得到了这样一个分布的均值(上面的紫色图表/黄色星星)。问题是是否存在一个合理的公式将正态分布的方差转化为 Beta 分布方差的等效形式。为了回答这个问题,我随机生成了几千个具有已知均值和方差的正态技能差异分布,然后将这些随机生成的分布中的每个点通过 sigmoid 函数转换为其 0 到 1 之间的等效概率。以下是一些示例,其中我调整了上面紫色分布的方差,保持相同的均值:
随着方差的增加,已知的信息减少,我们对蓝色玩家比红色玩家更优秀的信心也减少。黄色是方差最高的分布。
我在 Beta 旁边加了一个问号,因为没有保证通过这样做生成 Beta 分布。实际上,我确实没有。然而,仅仅通过目测这些分布,我猜测我可以找到相应的 Beta 参数来很好地描述这些分布。让我们使用最大似然估计经验性地找到上述分布的参数,以拟合最佳的 Beta 参数。
最小化负对数似然,我得出了以下结果:
Beta 参数拟合远非完美,但在方差较低或概率不极端的情况下似乎表现非常好。在生成了数千个这些后,我简单地使用PySR 包拟合了一个符号回归模型。我给它 Beta 均值和正态方差,看看它是否能预测 Beta 样本量(对应于方差)。然后我将得到 Beta 分布模型的所有参数。拟合模型后,PySR 给出了迄今为止最高的评分(基于简洁性和性能)的方程:
其中 mu_b 是 Beta 均值,sigma 平方是正态方差
我对这个方程的简单程度感到惊讶,对于这样一个简单的方程,结果居然如此好。虽然它比上面更远离完美,但你可以看到它对数据的拟合相当好:
这就是符号回归的美妙之处。我们不需要花费数分钟或数小时来优化成千上万的数据点,而是编写一个一行的函数,可以立即给出类似的结果。
将这个关键方程转化为实际的评级系统需要几个额外的步骤,这些步骤不是文章的重点。如果你感兴趣,这涉及使用 Beta 分布的贝叶斯更新,然后使用代数返回到正态分布。我还通过实验发现,由于比赛是成对比较(我们从两个玩家那里获取信息),实际中我们的有效样本量远远高于更新典型的 Beta 分布(如教科书中的掷硬币示例)。对于较高的样本量的调整只涉及添加一个平方项。以下是我用于实现它的所有代码:
当然,此时我仍然不确定我的实现有多好。它会比 TrueSkill 和 Glicko 差得多吗?我很可能只是浪费了大量时间。所以,我对 Alice 进行了之前相同的收敛测试:
(这次我花时间拟合 TrueSkill)
令我惊讶的是,它表现非常好。由于只需要少量基本操作,它可以每秒计算数千个评分,因此具有极高的可扩展性。该系统还具有无需大量参数调整的优点,同时保持了 Elo 系统的多样性和简洁性。当然,我想进一步推进,看看它在哪里崩溃:
当概率接近 0 或 1 时,预测结果与最大似然估计不匹配。因此,最终由于 Beta 分布根本上不是一个完美的拟合,这里存在一些劣势。我相信有一个更复杂的方程可以帮助解决这个问题,或者某种修正方法。也许一位统计学家会告诉我,有比 Beta 分布更适合建模这些结果的分布。然而,我对我创建的系统感到非常惊喜。没有任何系统能做到完美,总是有需要考虑的缺陷。当然,有几种方法可以将这些想法融入实际应用中——如果你感兴趣,可以关注我,我会继续写相关内容。我希望这样的成功能激励读者尝试符号回归以适应自己的使用案例。这里有一个链接到我创建所有这些可视化的笔记本。这是一个库的一部分,我希望在其中实现所有这些评分系统的向量化版本(尽管这可能还需要几个月时间)。
参考文献:
[1] Miles Cranmer, 用于科学的可解释机器学习:PySR 和 SymbolicRegression.jl (2023)
[2] Mark Glickman, 国际象棋等级评估综合指南 (1995)
[3] Ralf Herbrich, Tom Minka, Thore Graepel, TrueSkill™: 贝叶斯技能评分系统 (2007),神经信息处理系统第 20 卷
除非另有说明,否则所有图像均由作者提供。
使用 SQL 中的 HAVING 和 DISTINCT 子句
原文:
towardsdatascience.com/using-the-having-and-distinct-clauses-in-sql-d9e3be67b4be
你应该知道的两个重要 SQL 子句
·发表于 Towards Data Science ·4 分钟阅读·2023 年 1 月 25 日
–
SQL 是从数据库中提取数据的强大工具——无论是从一个表还是多个表中。
话虽如此,在有效分析数据时,有些子句特别重要。
本文讨论的两个子句是 HAVING 和 DISTINCT 子句。
为什么我们需要 HAVING 子句?
HAVING 子句的目的是在使用 GROUP BY 函数时充当 WHERE 子句的等效物。
如果你经常使用 SQL,你会知道 GROUP BY 子句对于在表中聚合值是非常重要的,例如,获得特定数据组的平均值,计算特定组的最大值或最小值——以及众多其他功能。
假设数据库中存在以下表格:
来源:作者使用 PostgreSQL 创建的表。表格显示在 pgAdmin4 中。
我们可以看到表格包含:
-
由字母表示的品牌
-
每个品牌的产品由数字 ID 表示
-
每种产品的价格
让我们假设我们希望确定表中每个品牌的平均产品价格,但仅在每个品牌在表中存在多个条目的情况下,即我们可以看到品牌 B 和 D 只有一个条目。因此,我们不希望这些条目包含在分析中。
我们如何使用 HAVING 子句计算在表中有多个条目的品牌的平均价格?
以下是该子句:
select brand, avg(price) from brand_table group by brand having count(brand)>1 order by brand;
执行此子句时,SQL 返回如下结果:
来源:表格由作者使用 PostgreSQL 创建。表格在 pgAdmin4 中显示。
正如我们所见,SQL 只返回了品牌 A、C 和 E 的平均价格。由于品牌 B 和 D 在表中只有一个条目,因此这些品牌没有被包括在内。
如果我们仅选择使用 GROUP BY 语句而不使用 HAVING 子句,SQL 将包括所有品牌:
来源:表格由作者使用 PostgreSQL 创建。表格在 pgAdmin4 中显示。
然而,B 和 D 的上述价格并不特别重要——因为计算只有一个产品的平均价格没有意义。因此,我们使用 HAVING 子句来仅显示多个条目的品牌。
DISTINCT 子句的目的
最简单来说,DISTINCT 子句的目的是选择表中唯一的条目而没有重复。
在上表中,我们可以看到每个品牌条目有多个产品。然而,如果我们只是想在表中显示每个品牌,我们可以如下使用 DISTINCT:
select distinct(brand) from brand_table order by brand;
来源:数据由作者使用 PostgreSQL 生成。表格在 pgAdmin 4 中显示。
在上表中,你会注意到品牌 E 的两个产品价格都是 10。
当选择品牌和价格而不使用 DISTINCT 子句时,我们可以看到这两个条目都被显示:
select brand, price from brand_table order by brand, price;
来源:数据由作者使用 PostgreSQL 生成。表格在 pgAdmin 4 中显示。
然而,当我们包含 DISTINCT 子句时——我们看到表格仅显示了品牌 E 的一个条目,其中价格为 10——这正是我们想要的:
来源:数据由作者使用 PostgreSQL 生成。表格在 pgAdmin 4 中显示。
从这个角度来看,使用 DISTINCT 可以在必要时获取特定组中的唯一值——在这个案例中是品牌。即使价格相同的产品 ID 可能是唯一的——但在按品牌分析时这并不相关。
结论
在本文中,你已经看到了:
-
如何在 GROUP BY 语句中将 HAVING 用作 WHERE 等效项
-
使用 DISTINCT 来返回特定类别中的唯一值
非常感谢阅读,如有任何问题或反馈,我们将不胜感激!
免责声明:本文是“按原样”编写的,不提供任何担保。它的目的是提供数据科学概念的概述,不应被解读为专业建议。本文中的发现和解释仅代表作者的观点,不代表或与本文中提到的任何第三方相关联。作者与本文中提到的任何第三方没有关系。
在 Python 中使用 Tqdm 与 Asyncio
原文:
towardsdatascience.com/using-tqdm-with-asyncio-in-python-5c0f6e747d55
PYTHON CONCURRENCY
一种有效监控并发任务进度的方法
·发表于 Towards Data Science ·6 分钟阅读·2023 年 5 月 2 日
–
图片由 Jungwoo Hong 提供,来自 Unsplash
介绍
让我烦恼的事
对数据科学家来说,使用 Python 进行并发编程以提高效率并不罕见。观察后台各种子进程或并发线程,以保持计算或 IO 绑定任务的有序总是令人满意的。
但让我困扰的一件事是,当我在后台并发处理数百或数千个文件或执行数百个进程时,我总是担心是否有几个任务会悄悄挂起,从而使整个代码永远无法完成。我也很难知道代码目前的执行状态。
最糟糕的部分是,当我看着空白屏幕时,很难判断我的代码还需要多长时间才能执行完成或预计时间。这对我组织工作进度的能力非常有害。
因此,我想要一种方法来让我知道代码执行到了哪里。
过去是怎么做的
更传统的方法是在线程之间共享一个内存区域,在这个内存区域中放置一个计数器,让这个计数器在任务完成时加 1,然后使用一个线程不断打印这个计数器的值。
这不是一个好的解决方案:一方面,我需要在现有业务逻辑中添加计数代码,这违反了“低耦合,高内聚”的原则。另一方面,由于线程安全问题,我不得不非常小心锁机制,这会导致不必要的性能问题。
tqdm 就是解决方案
tqdm 使用进度条来指示任务的进度。图片由作者提供
有一天,我发现了 tqdm 库,它使用进度条来可视化我的代码进度。我可以用进度条来可视化我的 asyncio 任务的完成情况和预计完成时间吗?
我继续研究,最终成功了。然后我与大家分享这种方法,以便每个程序员都有机会监控他们的并发任务进度。让我们开始吧。
Python 中 asyncio 的背景
在我们开始之前,我希望你对 Python asyncio 有一些背景了解。我的文章描述了一些 asyncio 常见 API 的用法,这将帮助我们更好地理解 tqdm 的设计:
asyncio.gather、asyncio.as_completed 和 asyncio.wait 的最佳实践
tqdm 概览
正如官方网站所描述,tqdm 是一个用于显示循环进度条的工具。它使用简单,高度可定制,并且资源占用极低。
一种典型的用法是将一个可迭代对象传递给 tqdm 构造函数,然后你将得到一个如下所示的进度条:
或者你可以手动遍历并更新进度条的进度,当文件被读取时:
使用 tqdm 指示读取大数据集的进度。图片由作者提供
将 tqdm 与 asyncio 集成
总体来说,tqdm 使用起来非常简单。然而,GitHub 上关于将 tqdm 与 asyncio 集成的信息还不够多。因此,我深入挖掘了源代码,以查看 tqdm 是否支持 asyncio。
幸运的是,tqdm 的最新版本提供了 tqdm.asyncio
包,该包提供了 tqdm_asyncio
类。
tqdm_asyncio
类有两个相关的方法。一个是 tqdm_asyncio.as_completed
。从源代码中可以看出,它是 asyncio.as_completed
的一个封装:
@classmethod
def as_completed(cls, fs, *, loop=None, timeout=None, total=None, **tqdm_kwargs):
"""
Wrapper for `asyncio.as_completed`.
"""
if total is None:
total = len(fs)
kwargs = {}
if version_info[:2] < (3, 10):
kwargs['loop'] = loop
yield from cls(asyncio.as_completed(fs, timeout=timeout, **kwargs),
total=total, **tqdm_kwargs)
另一个是 tqdm_asyncio.gather
,从源代码中可以看出,它基于对 tqdm_asyncio.as_completed
的实现,该实现模拟了 asyncio.gather
的功能:
@classmethod
async def gather(cls, *fs, loop=None, timeout=None, total=None, **tqdm_kwargs):
"""
Wrapper for `asyncio.gather`.
"""
async def wrap_awaitable(i, f):
return i, await f
ifs = [wrap_awaitable(i, f) for i, f in enumerate(fs)]
res = [await f for f in cls.as_completed(ifs, loop=loop, timeout=timeout,
total=total, **tqdm_kwargs)]
return [i for _, i in sorted(res)]
接下来,我将描述这两个 API 的用法。在我们开始之前,我们还需要做一些准备工作。在这里,我写了一个简单的方法,用于模拟一个具有随机睡眠时间的并发任务:
紧接着,我们将创建 2000 个并发任务,然后使用 tqdm_asyncio.gather
替代熟悉的 asyncio.gather
方法,看看进度条是否正常工作:
tqdm_asyncio.gather
的效果。图片由作者提供
噠噠!我终于知道我的任务完成了。相当酷。
或者让我们用 tqdm_asyncio.as_completed
替换 tqdm_asyncio.gather
并再试一次:
tqdm_asyncio.as_completed 也运行良好。图像来源:作者
很好,它仍然运行良好。
高级技巧和窍门
一些常见的配置项
tqdm 有一套丰富的 配置项,这里列出一些常见的配置项。
desc
。你可以配置一个 desc 参数,在进度条前显示标题,这在区分多个任务组时很有用。
desc 配置项的作用。图像来源:作者
ncols
。如果默认的进度条过短,可以通过这个参数将其延长。
使用 ncols 更改进度条的宽度。图像来源:作者
colour
。PyCharm 的 CLI 默认将进度条显示为红色,这仍然有些刺眼,因此你可以使用这个参数将进度条更改为其他颜色。但截至写作本文时,我仍未找到将文本更改为白色的方法。
使用 color 来更改进度条的颜色。图像来源:作者
bar_format
。这个选项允许你灵活控制进度条显示的内容和格式。例如,如果你想在顶部显示 ETA。
使用 bar_format 自定义进度条的内容。图像来源:作者
异常处理
从源代码中可以看出,tqdm 通过 tqdm_asyncio.as_completed
方法实现了 gather
方法。因此,我们不能通过使用 return_exceptions
参数来跳过异常捕获。
这很遗憾。但我们仍然可以在 tqdm_asyncio.as_completed
中通过 try…exception
处理异常:
异常处理。图像来源:作者
现实世界中的应用案例
许多 asyncio 的代码示例使用 asyncio.sleep
来模拟 IO 绑定的情况,但这不幸地过于简化了现实情况。我们应该使用实际案例来解释如何在 asyncio 中使用 tqdm。
然而,由于篇幅原因,我们无法在本章中使用实际案例。在下一章中,我们将演示 tqdm 进度条如何在实际使用 asyncio 实现 map-reduce 程序处理大型文件的示例中工作。
## 在 Python 中结合多进程和 asyncio 提升性能
使用实际案例来解释代码实现
towardsdatascience.com
结论
在 asyncio 代码中使用 tqdm 指示进度有很多好处:
-
我们可以在调用者的进度条中显示进度,而无需干涉业务代码。
-
所有工作都可以在主进程中完成,无需担心线程安全和性能问题。
-
图形化展示总是比枯燥的文字描述要生动得多。
-
所有这些只需一行代码。
我还尝试过其他进度条库,比如 alive-progress,它在展示效果上要酷得多,但 alive-progress 不支持 asyncio。
如果设置得当,tqdm 也可以产生一些很酷的效果,但由于时间原因我没有深入研究,所以欢迎进一步讨论并留下评论。你可能会帮助到更多感兴趣的读者。
通过 加入 Medium,你将可以无限制地访问我所有的帖子以及成千上万其他作者的文章。这只需花费你一杯咖啡的钱,但对我来说是极大的鼓励。
本文最初发布在: www.dataleadsfuture.com/using-tqdm-with-asyncio-in-python/
使用 pykrige 和 matplotlib 进行地质变化的空间可视化
探索井测量数据中的空间地质变化
· 发布于 Towards Data Science ·阅读时长 7 分钟·2023 年 6 月 13 日
–
挪威大陆架上声波压缩缓慢度测量的空间变化。图像由作者提供。
在处理地质和岩石物理数据时,我们通常希望了解这些数据在我们研究的区域或场地上的变化。我们可以通过网格化实际测量值,并推测尚未通过钻孔探索的其他区域的数据来实现这一点。
一种进行这种外推的方法是克里金法,这是一种以南非矿业工程师 Danie G. Krige 命名的地统计学程序。克里金法的核心思想在于其估计技术:利用观察数据之间的空间相关性来预测未测量位置的值。
通过衡量变量在距离上的变化,这种方法建立了一种统计关系,可以用于预测区域内的值,将分散的数据点转化为连贯的空间地图。
在本教程中,我们将探讨一个名为pykrige的 Python 库。该库设计用于二维和三维克里金计算,并且使用井测数据非常简便。
导入库和数据
首先,我们需要导入我们将要用到的库。对于这篇文章,我们将需要以下库:
-
pandas — 用于读取我们以
csv
格式存储的数据 -
matplotlib 用于创建我们的可视化
-
pykrige 用于执行克里金计算
-
numpy 用于一些数值计算
import pandas as pd
import matplotlib.pyplot as plt
from pykrige import OrdinaryKriging
import numpy as np
一旦我们导入了库,现在我们可以导入我们的数据。
在本教程中,我们将使用从 Xeek 和 Force 2020 机器学习竞赛中获得的用于预测井日志测量的岩性数据集。该数据集的详细信息可以在本文底部找到。
这个竞赛数据集的子集包含 65 个井位置及 Balder 组的平均声学压缩慢度测量值。
要读取数据,我们可以使用 pandas 的 read_csv()
函数,并传入数据文件的位置。在此示例中,我们使用相对于 Jupyter Notebook 的路径,但如果文件位于其他位置,我们也可以使用绝对路径。
df = pd.read_csv('Data/Xeek Force 2020/Xeek_2020_Balder_DTC_AVG.csv')
df
当我们查看数据框时,我们将看到我们有 65 口井,包含 Balder 组顶部的位置(X_LOC 和 Y_LOC 为网格坐标,LAT 和 LON 为纬度和经度)。我们还记录了遇到该地层的真实垂直深度(TVDSS)以及声学压缩慢度的均值(DTC)。
包含我们选择的井位置数据和挪威北海 Balder 组的 DTC — 声学压缩慢度值的数据框。图片来源于作者。
可视化井的空间位置
现在我们的数据已成功加载到数据框中,我们可以可视化数据,以了解井的位置。为此,我们将使用 matplotlib 的散点图,并传入经度和纬度列。
plt.scatter(df['Longitude'], df['Latitude'], c=df['DTC'])
当我们运行上述代码时,我们会得到以下图表。
基本的 matplotlib 图形展示了我们在挪威北海区域的井位置和 DTC 值。图片来源于作者。
我们可以看到上面的图形非常基本,没有颜色条或坐标轴标签。
让我们稍微修改图表,添加这些特性。
cm = plt.cm.get_cmap('viridis')
plt.figure(figsize=(10,10))
scatter = plt.scatter(df['LON'], df['LAT'], c=df['DTC_MEAN'], cmap=cm, s=50)
plt.colorbar(scatter)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.show()
当我们运行上述代码时,我们会得到以下图形,这为我们提供了更多有关数据的信息。我们可以使用颜色条来估算点值。
在添加颜色条和标签后的挪威北海选择井的 matplotlib 散点图。图片来源于作者。
应用克里金法
为了更好地理解数据点及 DTC 测量在 Balder 组区域的变化,我们可以使用克里金法和数据点来填补测量值之间的空白。
为此,我们需要从 pykrige 库中创建一个 OrdinaryKriging
对象。
在此对象中,我们传入了 x 和 y 的位置数据,以及我们要映射到 z 参数的数据。
我们还需要选择要使用的变差函数模型。在此案例中,我们将使用指数模型。有关模型类型的更多细节,请参见文档。
由于我们使用纬度和经度作为 x 和 y 坐标,我们可以将 coordinates_type 参数更改为 geographic
OK = OrdinaryKriging(x=df['LON'],
y=df['LAT'],
z=df['DTC_MEAN'],
variogram_model='exponential',
verbose=True, enable_plotting=True,
coordinates_type='geographic')
当我们运行上述代码时,我们会返回以下模型总结和半变差函数。
来自 pykrige 的模型总结。图片由作者提供。
下面是返回的参数的简要总结:
-
颗粒:颗粒是变差函数的 y 截距,表示零距离处的方差,通常是由于测量误差或非常小尺度的变化造成的。
-
完全天花板:天花板是变差函数达到的最大方差并开始趋于平稳的点,这发生在点之间距离非常远时。
-
范围:范围是变差函数达到天花板的距离,意味着在此距离之外,进一步分离点不会增加方差。
-
部分天花板:部分天花板是天花板和颗粒之间的差异,表示数据中空间结构的方差量。
这可以让我们理解模型根据生成的线条和点的形状对数据的适用性。
显示克里金结果
要开始显示我们的数据,我们需要创建一个数据网格。
为此,我们首先创建纬度和经度的数组,介于我们定义的坐标之间。在这种情况下,我们希望图表从北纬 57.5 度扩展到北纬 62 度,从东经 1.5 度扩展到东经 4.5 度。
使用 np.arange
将允许我们以规则间隔创建这些数组。
grid_lat = np.arange(57.5, 62, 0.01, dtype='float64')
grid_long = np.arange(1.5, 4.5, 0.01,dtype='float64')
现在我们有了 X 和 Y 坐标,我们可以创建我们的值网格。为此,我们调用 OK.execute
,并传入我们的纬度和经度数组。
zstar, ss = OK.execute('grid', grid_long, grid_lat)
这将返回两个数组。我们的数据网格(zstar)和与之相关的不确定性(ss)
接下来,我们可以使用我们的数据数组,并使用 matplotlib 的 imshow
绘制它。
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10,10))
image = ax.imshow(zstar, extent=(1.5, 4.5, 57.5, 62), origin='lower')
ax.set_xlabel('Longitude', fontsize=14, fontweight='bold')
ax.set_ylabel('Latitude', fontsize=14, fontweight='bold')
scatter = ax.scatter(x=df['LON'], y=df['LAT'], color='black')
colorbar = fig.colorbar(image)
colorbar.set_label('DTC (us/ft)', fontsize=14, fontweight='bold')
plt.show()
当我们运行此代码时,我们得到以下地图,显示了在我们 65 口井中 Balder 地层的声学压缩慢度的变化。
使用 pykrige 生成的声学压缩慢度 (DTC) 数据网格。图片由作者提供。
我们可以看到,在北纬 59 至 60 度之间,我们有更快的岩石,而在东北和西南地区,我们有更慢的岩石。
为了解释这一点,我们需要了解每个井的位置的形成深度。这将帮助我们确定差异是否与埋藏和压实或其他地质过程有关。
我们将在未来的文章中看到如何做到这一点。
可视化克里金不确定性
在查看此类数据时,一个关键点是理解与克里金相关的不确定性。
我们可以通过重复使用相同的绘图代码来实现这一点,而不是传入 zstar
,我们可以将其替换为之前创建的 ss
变量。
fig, ax = plt.subplots(figsize=(10,10))
image = ax.imshow(ss, extent=(1.5, 4.5, 57.5, 62), origin='lower')
ax.set_xlabel('Longitude', fontsize=14, fontweight='bold')
ax.set_ylabel('Latitude', fontsize=14, fontweight='bold')
scatter = ax.scatter(x=df['LON'], y=df['LAT'], color='black')
colorbar = fig.colorbar(image)
colorbar.set_label('DTC (us/ft)', fontsize=14, fontweight='bold')
plt.show()
通过以下图表,我们能够看到不确定性高或低的区域。
由 pykrige 生成的声波压缩慢度(DTC)不确定性数据网格。图像由作者提供。
在我们井覆盖较少的区域,我们的不确定性会更高,而在我们有多个井的区域,我们的不确定性会更低。
总结
在本教程中,我们已经看到如何对井日志测量值(DTC)进行平均处理,并将其映射到整个区域。这使我们能够了解地理区域内数据的趋势。
然而,在查看这些数据时,我们必须记住,我们看到的是一个 2D 表面,而不是我们在地下遇到的更复杂的 3D 结构。因此,测量的变化可能是由深度变化造成的。
使用的数据集
本文使用的数据集是由 Xeek 和 FORCE 2020 (Bormann 等,2020) 组织的机器学习比赛使用的训练数据集的一个子集。它在挪威政府发布的 NOLD 2.0 许可证下发布,详细信息可以在这里找到:挪威开放政府数据许可证(NLOD)2.0。完整数据集可以在这里访问。
数据集的完整参考资料是:
Bormann, Peter, Aursand, Peder, Dilib, Fahad, Manral, Surrender, & Dischington, Peter. (2020). FORCE 2020 Well well log and lithofacies dataset for machine learning competition [数据集]. Zenodo. doi.org/10.5281/zenodo.4351156
感谢阅读。在你离开之前,你应该一定要订阅我的内容,并将我的文章发送到你的邮箱中。 你可以在这里完成订阅!
其次,你可以通过订阅会员来获得完整的 Medium 体验,并支持成千上万的其他作家和我。它每月只需$5,并且你可以完全访问所有精彩的 Medium 文章,还可以有机会通过你的写作赚钱。
如果你通过 我的链接注册,你将通过你的一部分费用直接支持我,这不会增加你的费用。如果你这样做,非常感谢你的支持。
利用 PyArrow 改进 pandas 和 Dask 工作流
现在最大限度地利用 PyArrow 在 pandas 和 Dask 中的支持
·
关注 发表在 Towards Data Science · 13 min 阅读 · 2023 年 6 月 6 日
–
所有图片均由作者创建
介绍
本文探讨了我们现在可以在何处使用 PyArrow 来改进我们的 pandas 和 Dask 工作流。对 PyArrow 数据类型的总体支持在 pandas 2.0 中添加到 pandas 和 Dask 中。这解决了两个库用户长期以来的一些痛点。pandas 用户常常向我抱怨 pandas 不支持任意数据类型中的缺失值,或非标准数据类型支持不够好。Dask 用户特别烦恼的问题是处理大数据集时内存不足。与 NumPy 对象列相比,PyArrow 支持的字符串列可以减少多达 70% 的内存消耗,因此有潜力缓解这个问题,并提供巨大的性能提升。
pandas 和 Dask 对 PyArrow 数据类型的支持仍然相对较新。建议在至少 pandas 2.1 发布之前,谨慎选择 PyArrow 的 dtype_backend
。并非所有 API 部分都已优化。然而,在某些工作流中,你应该能获得显著的改进。本文将介绍几个我建议立即切换到 PyArrow 的示例,因为它已经提供了巨大的好处。
Dask 本身可以通过 PyArrow 数据类型在多方面受益。我们将探讨 PyArrow 支持的字符串如何轻松缓解 Dask 集群内存不足的问题,以及我们如何通过利用 PyArrow 改善性能。
我是 pandas 核心团队的一员,并且在实现和改进 pandas 中的 PyArrow 支持方面参与了大量工作。我最近加入了 Coiled,在那里我负责 Dask。我的任务之一是改进 Dask 中的 PyArrow 集成。
PyArrow 支持的总体概述
PyArrow 数据类型最初是在 pandas 1.5 中引入的。该实现是实验性的,我不建议在 pandas 1.5.x 上使用它。对这些数据类型的支持仍然相对较新。pandas 2.0 提供了巨大的改进,包括使选择 PyArrow 支持的 DataFrames 变得容易。我们仍在努力在各处提供适当的支持,因此在至少 pandas 2.1 发布之前应谨慎使用。这两个项目都在不断努力改进对 Dask 和 pandas 的支持。
我们鼓励用户尝试使用它们!这将帮助我们更好地了解仍然缺乏支持或不够快的地方。反馈有助于我们改进支持,并将显著减少创建流畅用户体验所需的时间。
数据集
我们将使用包含所有 Uber 和 Lyft 乘车记录的纽约市出租车数据集。它具有一些有趣的属性,如价格、小费、司机薪酬等。数据集可以在这里找到(参见 服务条款),并存储在 parquet 文件中。在分析 Dask 查询时,我们将使用一个公开的 S3 存储桶以简化我们的查询:s3://coiled-datasets/uber-lyft-tlc/
。我们将使用 2022 年 12 月的数据集进行 pandas 查询,因为这是在我的机器(24GB 内存)中舒适地适配的最大数据集。我们必须避免过度使用 RAM,因为这可能会在分析性能时引入副作用。
我们还将研究 read_csv
的性能。我们将使用可以在这里找到的 芝加哥犯罪 数据集。
Dask 集群
设置 Dask 集群有各种不同的选项,参见 Dask 文档 获取非详尽的部署选项列表。我将使用 Coiled 在 AWS 上通过 30 台机器创建一个集群。
import coiled
cluster = coiled.Cluster(
n_workers=30,
name="dask-performance-comparisons",
region="us-east-2", # this is the region of our dataset
worker_vm_type="m6i.large",
)
Coiled 已连接到我的 AWS 账户。它在我的账户中创建集群并为我管理所有资源。30 台机器足以舒适地操作我们的数据集。我们将研究如何通过一些小的修改将所需的工作节点数减少到 15 台。
pandas StringDtype
由 PyArrow 支持
我们从 pandas 1.0 中最初引入的一个功能开始。将 pandas 或 Dask 中的 dtype 设置为 string
会返回一个具有 StringDtype
的对象。这个功能相对成熟,应该能提供平滑的用户体验。
历史上,pandas 通过 dtype 为 object
的 NumPy 数组来表示字符串数据。NumPy 对象数据存储为指向内存中实际数据的指针数组。这使得遍历包含字符串的数组非常缓慢。pandas 1.0 最初引入了所谓的 StringDtype
,允许对字符串进行更简单和一致的操作。这个 dtype 仍然由 Python 字符串支持,因此性能也不佳。它提供了字符串数据的清晰抽象。
pandas 1.3 最终引入了一种创建高效字符串 dtype 的增强功能。该数据类型由 PyArrow 数组支持。PyArrow 提供了一种数据结构,能够进行高效且节省内存的字符串操作。从那时起,用户可以使用在内存中连续的字符串 dtype,因此非常快速。该 dtype 可以通过 string[pyarrow]
请求。或者,我们可以通过指定 string
作为 dtype 并设置来请求它:
pd.options.mode.string_storage = "pyarrow"
由于 Dask 构建在 pandas 之上,这里也可以使用这种字符串数据类型。除此之外,Dask 还提供了一个方便的选项,可以自动将所有字符串数据转换为 string[pyarrow]
。
dask.config.set({"dataframe.convert-string": True})
这是避免字符串列使用 NumPy 对象数据类型的方便方法。此外,它还具有创建原生 PyArrow 数组的优势,适用于处理 Arrow 对象的 I/O 方法。在提供显著性能提升的同时,PyArrow 字符串消耗的内存显著减少。与 NumPy 对象相比,使用 PyArrow 字符串的平均 Dask DataFrame 的内存消耗约为原始内存的 33-50%。这解决了 Dask 用户在处理大型数据集时内存不足的最大痛点。该选项启用了 Dask 测试套件中的全局测试。这确保了 PyArrow 支持的字符串足够成熟,以提供流畅的用户体验。
让我们看几个典型的字符串操作。我们将从一些 pandas 示例开始,然后切换到在 Dask 集群上执行的操作。
我们将使用 df.convert_dtypes
将我们的对象列转换为 PyArrow 字符串数组。还有更高效的方法来获取 pandas 中的 PyArrow 数据类型,我们会在后面探讨。我们将使用 2022 年 12 月的 Uber-Lyft 数据集,这个文件在我的机器上可以舒适地放入内存中。
import pandas as pd
pd.options.mode.string_storage = "pyarrow"
df = pd.read_parquet(
"fhvhv_tripdata_2022-10.parquet",
columns=[
"tips",
"hvfhs_license_num",
"driver_pay",
"base_passenger_fare",
"dispatching_base_num",
],
)
df = df.convert_dtypes(
convert_boolean=False,
convert_floating=False,
convert_integer=False,
)
在这个示例中,我们的 DataFrame 对所有非字符串列都有 NumPy 数据类型。让我们开始筛选所有由 Uber 操作的乘车记录。
df[df["hvfhs_license_num"] == "HV0003"]
这个操作创建了一个 True/False 值的掩码,指定 Uber 是否进行了乘车。这不使用任何特殊的字符串方法,但等式比较会交给 PyArrow。接下来,我们将使用 pandas 实现的 String 访问器,它允许你对每个元素执行各种字符串操作。我们想要找到所有从以 "B028"
开头的基地发出的乘车记录。
df[df["dispatching_base_num"].str.startswith("B028")]
startswith
会遍历我们的数组,检查每个字符串是否以指定的子字符串开头。PyArrow 的优势显而易见。数据在内存中是连续的,这意味着我们可以高效地进行遍历。此外,这些数组还具有一个第二个数组,其中包含指向每个字符串首个内存地址的指针,这使得计算起始序列更快。
最后,我们查看一个 GroupBy
操作,它对 PyArrow 字符串列进行分组。分组的计算也可以交给 PyArrow,这比在 NumPy 对象数组上因子化要高效得多。
df.groupby(
["hvfhs_license_num", "dispatching_base_num"]
).mean(numeric_only=True)
让我们看看这些操作与字符串列由 NumPy 对象数据类型表示的 DataFrames 的对比。
结果或多或少符合我们的预期。基于字符串的比较在 PyArrow 字符串上执行时显著更快。大多数字符串访问器应该能提供巨大的性能提升。另一个有趣的观察是内存使用,与 NumPy 对象 dtype 相比,减少了大约 50%。我们将进一步用 Dask 详细研究这一点。
Dask 镜像了 pandas API,并将大多数操作委派给 pandas。因此,我们可以使用相同的 API 来访问 PyArrow 字符串。上述选项是全局请求这些选项的便利方式,我们将在这里使用:
dask.config.set({"dataframe.convert-string": True})
在开发过程中,这种选项的最大好处之一是能够轻松地在 Dask 中全局测试 PyArrow 字符串,以确保一切顺利。我们将利用 Uber-Lyft 数据集进行探索。该数据集在我们的集群上占用大约 240GB 的内存。我们的初始集群有 30 台机器,足以舒适地进行计算。
import dask
import dask.dataframe as dd
from distributed import wait
dask.config.set({"dataframe.convert-string": True})
df = dd.read_parquet(
"s3://coiled-datasets/uber-lyft-tlc/",
storage_options={"anon": True},
)
df = df.persist()
wait(df) # Wait till the computation is finished
我们将数据保存在内存中,以避免 I/O 性能对我们性能测量的影响。我们的数据现在已存储在内存中,这使得访问速度很快。我们将执行类似于 pandas 计算的操作。主要目标之一是展示 pandas 的好处如何转化为在 Dask 分布式环境中的计算。
第一个观察结果是,使用 PyArrow 支持的字符串列的数据框仅消耗 130GB 内存,仅为使用 NumPy 对象列时的一半。我们的数据框中只有几个字符串列,这意味着在切换到 PyArrow 字符串时,字符串列的内存节省实际上比约 50% 更高。因此,我们将在对 PyArrow 字符串列执行操作时将集群大小减少到 15 个工作节点。
cluster.scale(15)
我们通过对数据框进行后续筛选来衡量掩码操作和一个字符串访问器的性能。
df = df[df["hvfhs_license_num"] == "HV0003"]
df = df[df["dispatching_base_num"].str.startswith("B028")]
df = df.persist()
wait(df)
我们可以看到,可以使用与之前示例中相同的方法。这使得从 pandas 迁移到 Dask 相对容易。
此外,我们将再次对数据执行 GroupBy
操作。这在分布式环境中要困难得多,因此结果更具趣味性。之前的操作在大型集群上相对容易并行,而 GroupBy
则更具挑战性。
df = df.groupby(
["hvfhs_license_num", "dispatching_base_num"]
).mean(numeric_only=True)
df = df.persist()
wait(df)
我们得到了 2 和 3 倍的显著改进。这一点尤其引人注目,因为我们将集群的大小从 30 台机器减少到 15 台,成本降低了 50%。随后,我们还将计算资源减少了 2 倍,这使得我们的性能提升更加令人印象深刻。因此,性能分别提高了 4 倍和 6 倍。我们可以在较小的集群上执行相同的计算,这样既节省了成本,也在总体上更高效,并且仍然能够获得性能提升。
总结一下,我们看到 PyArrow 字符串列在与 DataFrame 中的 NumPy 对象列比较时有了巨大的改进。切换到 PyArrow 字符串是一个相对较小的改变,但可能极大地提升依赖于字符串数据的平均工作流程的性能和效率。这些改进在 pandas 和 Dask 中都同样明显!
I/O 方法中的 Engine 关键字
我们现在将看看 pandas 和 Dask 中的 I/O 函数。一些函数有自定义实现,如 read_csv
,而其他函数则委派给另一个库,如 read_excel
委派给 openpyxl
。其中一些函数新增了 engine
关键字,使我们能够委派给 PyArrow
。PyArrow 解析器默认是多线程的,因此可以提供显著的性能提升。
pd.read_csv("Crimes_-_2001_to_Present.csv", engine="pyarrow")
这个配置将返回与其他引擎相同的结果。唯一的区别是使用了 PyArrow 来读取数据。read_json
也提供了相同的选项。添加 PyArrow-engines 是为了提供更快的数据读取方式。改进的速度只是一个优势。PyArrow 解析器返回的数据是一个 PyArrow Table。PyArrow Table 提供了内置功能来转换为 pandas DataFrame
。根据数据的不同,这可能需要在转换为 NumPy(字符串、带缺失值的整数等)时进行拷贝,从而带来不必要的减速。这就是 PyArrow dtype_backend
的作用。它在 pandas 中实现为一个 ArrowExtensionArray
类,背后是一个 PyArrow ChunkedArray。直接的结果是,从 PyArrow Table 转换为 pandas 非常便宜,因为它不需要任何拷贝。
pd.read_csv(
"Crimes_-_2001_to_Present.csv",
engine="pyarrow",
dtype_backend="pyarrow",
)
这将返回一个由 PyArrow 数组支持的 DataFrame
。pandas 并未在所有方面都得到优化,因此这可能会在后续操作中造成减速。如果工作负载特别重 I/O,这可能是值得的。让我们看一下直接比较:
我们可以看到,与 C-engine 相比,PyArrow-engine 和 PyArrow dtypes 提供了 15 倍的加速。
相同的优势适用于 Dask。Dask 封装了 pandas 的 csv 读取器,因此,它免费获得了相同的功能。
对 Dask 的比较要复杂一些。首先,我的示例从本地机器读取数据,而我们的 Dask 示例将从 S3 存储桶中读取数据。网络速度将是一个相关因素。此外,分布式计算也有一些开销,我们需要考虑到这一点。
我们这里只关注速度,所以我们将从一个公共 S3 存储桶中读取一些时间序列数据。
import dask.dataframe as dd
from distributed import wait
df = dd.read_csv(
"s3://coiled-datasets/timeseries/20-years/csv/",
storage_options={"anon": True},
engine="pyarrow",
parse_dates=["timestamp"],
)
df = df.persist()
wait(df)
我们将为 engine="c"
、engine="pyarrow"
以及额外的 engine="pyarrow"
配置 dtype_backend="pyarrow"
执行这个代码片段。让我们来看看一些性能比较。两个示例都在集群上的 30 台机器上执行。
PyArrow 引擎运行速度约为 C 引擎的两倍。两种实现使用了相同数量的机器。使用 PyArrow dtype_backend
时内存使用量减少了 50%。如果仅将对象列转换为 PyArrow 字符串,也可以实现相同的减少,这能在后续操作中提供更好的体验。
我们已经看到 Arrow 引擎提供了显著的速度提升,超过了自定义 C 实现。它们尚未支持自定义实现的所有功能,但如果你的使用场景与支持的选项兼容,你应该能免费获得显著的速度提升。
使用 PyArrow dtype_backend
的情况稍微复杂一些。并非所有 API 区域都经过优化。如果你在 I/O 函数之外花费大量时间处理数据,那么这可能无法满足你的需求。如果你的工作流程在读取数据时花费大量时间,它会加速你的处理。
PyArrow 原生 I/O 读取器中的 dtype_backend
其他一些 I/O 方法也有一个 engine 关键字。read_parquet
是最常见的例子。不过这里的情况稍有不同。这些 I/O 方法默认已经使用 PyArrow 引擎。因此解析尽可能高效。另一个潜在的性能提升是使用dtype_backend
关键字。通常,PyArrow 会返回一个 PyArrow 表,然后转换为 pandas DataFrame。PyArrow 的数据类型会转换为其 NumPy 等效类型。设置dtype_backend="pyarrow"
可以避免这种转换。这能显著提高性能并节省大量内存。
让我们看看 pandas 性能的比较。我们从 2022 年 12 月读取了 Uber-Lyft 出租车数据。
pd.read_parquet("fhvhv_tripdata_2022-10.parquet")
我们读取了带有和不带有dtype_backend="pyarrow"
的数据。
我们可以很容易地看到,转换所花费的时间最多,读取 Parquet 文件后完成的转换。如果避免转换为 NumPy 数据类型,函数运行速度提高了三倍。
Dask 为read_parquet
提供了一个专门的实现,相比于 pandas 实现,具有一些针对分布式工作负载的优势。共同点是这两个函数都调度到 PyArrow 来读取 parquet 文件。两者都共同点是数据在成功读取文件后会转换为 NumPy 数据类型。我们正在读取整个 Uber-Lyft 数据集,这在我们的集群上消耗了约 240GB 的内存。
import dask.dataframe as dd
from distributed import wait
df = dd.read_parquet(
"s3://coiled-datasets/uber-lyft-tlc/",
storage_options={"anon": True},
)
df = df.persist()
wait(df)
我们在 3 种不同的配置下读取了数据集。首先是使用默认的 NumPy 数据类型,然后打开 PyArrow 字符串选项:
dask.config.set({"dataframe.convert-string": True})
最后是dtype_backend="pyarrow"
。让我们看看这在性能方面意味着什么:
类似于我们的 pandas 示例,我们可以看到转换为 NumPy dtypes 会占用我们大量的运行时间。PyArrow dtypes 提供了显著的性能提升。两个 PyArrow 配置所使用的内存是 NumPy dtypes 的一半。
PyArrow 字符串比一般的 PyArrow dtype_backend
更成熟。根据我们得到的性能图表,当使用 PyArrow 字符串和 NumPy dtypes 对于所有其他 dtypes 时,我们得到的性能提升大致相同。如果某个工作流程在 PyArrow dtypes 上的表现还不够好,我建议仅启用 PyArrow 字符串。
结论
我们已经看到如何在 pandas 和 Dask 中利用 PyArrow。PyArrow 支持的字符串列有潜力以积极的方式影响大多数工作流程,并为 pandas 2.0 提供顺畅的用户体验。Dask 提供了一个方便的选项,可以在可能的情况下全局避免使用 NumPy 对象 dtype,这使得选择 PyArrow 支持的字符串变得更加容易。PyArrow 在其他可用领域也提供了巨大的加速。PyArrow 的 dtype_backend
仍然相当新,但有潜力显著减少 I/O 时间。探索它是否能解决性能瓶颈是非常值得的。许多工作正在进行中,以改进对一般 PyArrow dtypes 的支持,并有潜力在不久的将来加速平均工作流程。
pandas 目前有一个提议,从 pandas 3.0 开始默认将字符串推断为 PyArrow 支持的字符串。此外,这还包括更多领域,其中更多地依赖 PyArrow 是非常有意义的(例如:Decimals、结构化数据等)。您可以在这里阅读提议的详细信息。
感谢阅读。欢迎在评论中分享您对两个库中 PyArrow 支持的想法和反馈。我将会写后续的文章,专注于这个话题以及 pandas 一般内容。如果您喜欢阅读关于 pandas 和 Dask 的更多内容,可以关注我在 Medium 上的更新。
V-Net,U-Net 在图像分割中的“大哥”
原文:
towardsdatascience.com/v-net-u-nets-big-brother-in-image-segmentation-906e393968f7
欢迎阅读这本关于 V-Net 的指南,它是著名的 U-Net 的“亲戚”,专用于 3D 图像分割。你将对它了如指掌!
·发表在Towards Data Science ·阅读时间 8 分钟·2023 年 7 月 28 日
–
欢迎踏上探索深度学习架构的激动人心的旅程!你可能已经对U-Net有所了解,它是计算机视觉领域的一个重大突破,显著重塑了图像分割的格局。
今天,让我们将焦点转向 U-Net 的“大哥”——V-Net。
由研究人员 Fausto Milletari、Nassir Navab 和 Seyed-Ahmad Ahmadi 发表的论文“VNet: Fully Convolutional Neural Networks for Volumetric Medical Image Segmentation”介绍了一种突破性的 3D 图像分析方法。
这篇文章将带你了解这篇开创性的论文,阐明其独特的贡献和架构进展。无论你是经验丰富的数据科学家、正在成长的 AI 爱好者,还是对最新科技趋势感兴趣的人,这里都有你需要的内容!
关于 U-Net 的简要提醒
在深入 V-Net 的核心之前,让我们先欣赏一下它的架构灵感——U-Net。如果这是你第一次接触 U-Net,别担心,我有一个关于 U-Net 架构的快速简单教程来帮助你。它如此简洁,你最多五分钟就能掌握这个概念!
下面是对 U-Net 的简要回顾:
U Net 架构,见于U Net 文章
U-Net 以其对称结构而闻名,呈现“U”字形。这种架构由两个不同的路径组成:
-
收缩路径(左): 在这里,我们逐渐减少图像的分辨率,同时增加过滤器的数量。
-
扩展路径(右): 这个路径作为收缩路径的镜像。我们逐渐减少过滤器的数量,同时增加分辨率,直到与原始图像大小对齐。
U-Net 的美在于其创新使用的**“残差连接”或“跳跃连接”**。这些连接收缩路径和扩展路径中的对应层,使网络能够保留在收缩过程中通常丢失的高分辨率细节。
残差连接,来自 U Net 论文
为什么这很重要?因为它简化了反向传播过程中的梯度流动,特别是在早期层中。实质上,我们绕过了梯度消失的风险——一个常见的问题,梯度接近零,阻碍了学习过程:
图片来源:作者
现在,带着对 U-Net 的理解,让我们进入 V-Net 的世界。V-Net 的核心与 U-Net 相似,采用类似的编码器-解码器理念。但正如你很快会发现的,它具有一套独特的特征,使其与其兄弟 U-Net 有所区别。
V-Net 架构,来自 VNet 论文
V-Net 与 U-Net 有何不同?
让我们深入探讨吧!
区别 1:3D 卷积代替 2D 卷积
第一个区别显而易见。虽然 U-Net 是为 2D 图像分割量身定制的,但医学图像通常需要 3D 视角(例如体积脑扫描、CT 扫描等)。
这就是 V-Net 发挥作用的地方。V-Net 中的“V”代表“体积”,这种维度的变化要求将 2D 卷积替换为 3D 卷积。
区别 2:激活函数,PreLU 代替 ReLU
深度学习领域已经爱上了 ReLU 函数,因为它的简单性和计算效率。与 sigmoid 或 tanh 等其他函数相比,ReLU 是“非饱和”的,这意味着它减少了梯度消失的问题。
(左)ReLU,(中)LeakyReLU 和(右)PReLU,来自 PReLU 论文
但 ReLU 并不完美。它因一种被称为“Dying ReLU 问题”的现象而臭名昭著,其中许多神经元总是输出零,成为“死神经元”。为了解决这个问题,引入了 LeakyReLU,它在零的左侧有一个小但非零的斜率。
更进一步推理,V-Net 利用了参数化 ReLU(PReLU)。与其硬编码 LeakyReLU 的斜率,不如让网络来学习它?
毕竟,这是深度学习的核心哲学,我们希望尽可能少地引入归纳偏差,让模型自己学习一切,前提是我们有足够的数据。
区别 3:基于 Dice 分数的不同损失函数
现在,我们来到 V-Net 可能最具影响力的贡献——损失函数的改变。与 U-Net 的交叉熵损失函数不同,V-Net 使用 Dice 损失函数。
交叉熵函数,来自作者
但这个函数的主要问题是它对不平衡的类别处理不好。而这个问题在医学图像中非常常见,因为大多数时候背景的存在远多于感兴趣区域。
例如考虑这张图片:
背景无处不在,来自作者
结果是,一些模型可能变得“懒惰”,在任何地方都预测背景,因为它们仍然会得到一个较小的损失。
所以,V-Net 使用一种对这个问题更有效的损失函数:Dice 系数。
它更好的原因在于它将预测区域和真实值之间的重叠度量为比例,因此考虑了类别的大小。
尽管背景几乎无处不在,Dice 分数测量预测与真实值之间的重叠,因此即使类别占据主导地位,我们仍然可以得到一个介于 0 和 1 之间的数字。
Dice 系数,来自 VNet 论文
我认为这可能是本文的主要贡献,因为从 2D 卷积到 3D 卷积是处理 3D 图像的一个非常自然的想法。然而,这种损失函数已在图像分割任务中被广泛采用。
实际上,混合方法通常被证明是有效的,结合交叉熵损失和 Dice 损失,以利用两者的优点。
V-Net 的性能
因此,我们已经探讨了 V-Net 的独特方面,但你可能在想,“所有这些理论很棒,但 V-Net 实际上能否有效?”好吧,让我们对 V-Net 进行测试!
作者在 PROMISE12 数据集上评估了 V-Net 的性能。
PROMISE12 数据集是为 MICCAI 2012 前列腺分割挑战赛提供的。
V-Net 在 50 张磁共振(MR)图像上进行训练,这个数量并不多!
VNet 在 PROMISE 2012 挑战数据集上的分割,来自 VNet 论文
PROMISE 2012 挑战数据集上的定量指标,来自 VNet 论文
正如我们所见,即使标签很少,V-Net 也能够生成良好的定性分割,并获得非常好的 Dice 分数。
V-Net 的主要限制
确实,V-Net 在图像分割领域,特别是医学成像中设立了新的基准。然而,每项创新都有成长空间。在这里,我们将讨论 V-Net 可以改进的一些显著领域:
局限性 1:模型大小
从 2D 转到 3D 带来了显著的内存消耗增加。这种增加的连锁反应是多方面的:
-
该模型需要大量的内存空间。
-
它严重限制了批量大小(因为将多个 3D 张量加载到 GPU 内存中变得具有挑战性)。
-
医学成像数据稀缺且标注成本高,使得拥有如此多参数的模型更难以拟合。
局限性 2:不使用无监督学习或自监督学习
- V-Net 完全在监督学习的背景下运作,忽视了无监督学习的潜力。在未标注扫描大大超过标注扫描的领域中,融入无监督学习可能会带来突破性的改变。
局限性 3:没有不确定性估计
- V-Net 不估计不确定性,这意味着它无法评估自身预测的信心。这是贝叶斯深度学习闪耀的领域。(请参阅这篇贝叶斯深度学习简明介绍)。
局限性 4:缺乏鲁棒性
- 卷积神经网络(CNNs)传统上在泛化方面表现不佳。它们对于对比度变化、多模态分布或不同分辨率等变化不够鲁棒。这是 V-Net 还可以改进的另一个领域。
结论
V-Net,作为 U-Net 的较少知名但强大的对手,已经彻底改变了计算机视觉,尤其是图像分割。它从 2D 转向 3D 图像,并引入了现在广泛使用的 Dice 系数,为该领域设立了新的标准。
尽管存在局限性,V-Net 应该是开始进行 3D 图像分割任务的首选模型。为了进一步改进,探索无监督学习和整合注意力机制似乎是有前景的方向。
感谢阅读!在你离开之前:
- 查看我在 Github 上的AI 教程汇编
[## GitHub - FrancoisPorcher/awesome-ai-tutorials: 最佳 AI 教程汇编,助你成为…
最佳 AI 教程汇编,助你成为数据科学大师!- GitHub …
github.com](https://github.com/FrancoisPorcher/awesome-ai-tutorials?source=post_page-----906e393968f7--------------------------------)
你应该在你的收件箱中获取我的文章。 点击此处订阅。
如果你想访问 Medium 上的优质文章,只需每月$5 的会员订阅。如果你通过 我的链接进行注册,你可以在没有额外费用的情况下用你的一部分费用支持我。
如果你觉得这篇文章有见解且有帮助,请考虑关注我并点赞,以获取更多深入内容!你的支持帮助我继续制作有助于我们集体理解的内容。
参考文献
-
Milletari, F., Navab, N., & Ahmadi, S. A. (2016). V-Net:用于体积医学图像分割的全卷积神经网络。在 2016 年第四届国际 3D 视觉会议(3DV)上(第 565–571 页)。IEEE。
-
Ronneberger, O., Fischer, P., & Brox, T. (2015). U-Net:用于生物医学图像分割的卷积网络。在国际医学图像计算与计算机辅助干预会议(第 234–241 页)。Springer, Cham。