精通 NumPy 数值分析（二）
五、使用 NumPy 对批发分销商的客户进行聚类
通过查看 NumPy 在各种使用案例中的使用，您肯定可以提高自己的技能。 本章介绍的分析类型与迄今为止所见不同。 聚类是一种无监督的学习技术，用于理解和捕获数据集中的各种形式。 由于您没有标签来监督您的学习算法，因此在很多情况下，可视化是关键，这就是为什么您也会看到各种可视化技术的原因。
在本章中，我们将介绍以下主题：
- 无监督学习和聚类
- 超参数
- 扩展简单算法来聚类批发分销商的客户
无监督学习和聚类
让我们用一个例子快速回顾一下监督学习。 在训练机器学习算法时，您可以通过提供标签来观察和指导学习。 考虑下面的数据集，其中每一行代表一个客户，每一列代表一个不同的特征，例如年龄，性别，收入，职业 ，任期和城市。 看看这个表：
您可能需要执行不同类型的分析。 其中一项可能是预测哪些客户可能会离开，即客户流失分析。 为此，您需要根据每个客户的历史记录为其标记标签，以指示哪些客户已离开或留下，如下表所示：
您的算法将根据客户的标签了解客户的特征。 算法将学习离开或留下的客户的特征，并且，当您希望根据这些特征为新客户评分时，算法将根据该学习进行预测。 这称为监督学习。
关于此数据集可能要问的另一个问题可能是：此数据集中存在多少个客户组，以便其组中的每个客户都与其他客户相似，而与属于其他组的客户不同。 流行的聚类算法（例如 K 均值）可以帮助您解决这一问题。 例如，一旦 K 均值将客户分配到不同的集群，则一个集群可以大体上包括 30 岁以下且以 IT 作为职业的客户，而另一个集群可以大体上包括 60 岁以上的客户并且职业为老师。 您无需标记数据集即可执行此分析，因为算法足以查看记录并确定它们之间的相似性。 由于没有监督，因此这种学习称为无监督学习。
进行此类分析时，首先可视化数据集会很有帮助。 您可以从可用的数据集开始构建您的处理和建模工作流程。 以下代码段显示了如何使用plotly可视化三维数据集。plotly是一个库，可让您绘制许多不同的交互式图表以进行探索性分析，并使数据探索更加容易。
首先，您需要使用以下代码段安装必要的库：
## Installing necessary libraries with pip
!pip install plotly --user
!pip install cufflinks -user
然后，使用以下代码导入必要的库：
## Necessary imports
import os
import sys
import numpy as np
import pandas
import matplotlib.pyplot as plt
%matplotlib inline
import plotly.plotly as py
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
import cufflinks as cf
import plotly.graph_objs as go
init_notebook_mode(connected=True)
sys.path.append("".join([os.environ["HOME"]]))
您将使用sklearn.datasets模块中可用的iris数据集，如下所示：
from sklearn.datasets import load_iris
iris_data = load_iris()
iris数据具有四个特征； 它们如下：
iris_data.feature_names
['sepal length (cm)',
'sepal width (cm)',
'petal length (cm)',
'petal width (cm)']
首先，让我们检查一下前两个特征：
x = [v[0] for v in iris_data.data]
y = [v[1] for v in iris_data.data]
创建一个trace，然后创建数据和图形，如下所示：
trace = go.Scatter(
x = x,
y = y,
mode = 'markers'
)
layout= go.Layout(
title= 'Iris Dataset',
hovermode= 'closest',
xaxis= dict(
title= 'sepal length (cm)',
ticklen= 5,
zeroline= False,
gridwidth= 2,
),
yaxis=dict(
title= 'sepal width (cm)',
ticklen= 5,
gridwidth= 2,
),
showlegend= False
)
data = [trace]
fig= go.Figure(data=data, layout=layout)
plot(fig)
如下图所示，这将为您提供以下输出：
您可以继续查看其他变量，但是要更好地理解一张图表中要素之间的关系，可以使用scatterplot矩阵。 创建pandas.DataFrame与plotly结合使用会更加方便：
import pandas as pd
df = pd.DataFrame(iris_data.data,
columns=['sepal length (cm)',
'sepal width (cm)',
'petal length (cm)',
'petal width (cm)'])
df['class'] = [iris_data.target_names[i] for i in iris_data.target]
df.head()
使用plotly图形工厂，可以绘制scatterplot矩阵，如下所示：
import plotly.figure_factory as ff
fig = ff.create_scatterplotmatrix(df, index='class', diag='histogram', size=10, height=800, width=800)
plot(fig)
这将为您提供以下图表，如下图所示：
乍一看，花瓣长度，花瓣宽度和萼片长度似乎是建模的不错选择。 您可以通过以下代码使用 3D 图表进一步检查该数据集：
## Creating data for the plotly
trace1 = go.Scatter3d(
# Extracting data based on label
x=[x[0][0] for x in zip(iris_data.data, iris_data.target) if x[1] == 0],
y=[x[0][2] for x in zip(iris_data.data, iris_data.target) if x[1] == 0],
z=[x[0][3] for x in zip(iris_data.data, iris_data.target) if x[1] == 0],
mode='markers',
marker=dict(
size=12,
line=dict(
color='rgba(217, 217, 217, 0.14)',
width=0.5
),
opacity=0.8
)
)
trace2 = go.Scatter3d(
# Extracting data based on label
x=[x[0][0] for x in zip(iris_data.data, iris_data.target) if x[1] == 1],
y=[x[0][2] for x in zip(iris_data.data, iris_data.target) if x[1] == 1],
z=[x[0][3] for x in zip(iris_data.data, iris_data.target) if x[1] == 1],
mode='markers',
marker=dict(
color='rgb(#3742fa)',
size=12,
symbol='circle',
line=dict(
color='rgb(204, 204, 204)',
width=1
),
opacity=0.9
)
)
trace3 = go.Scatter3d(
# Extracting data based on label
x=[x[0][0] for x in zip(iris_data.data, iris_data.target) if x[1] == 2],
y=[x[0][2] for x in zip(iris_data.data, iris_data.target) if x[1] == 2],
z=[x[0][3] for x in zip(iris_data.data, iris_data.target) if x[1] == 2],
mode='markers',
marker=dict(
color='rgb(#ff4757)',
size=12,
symbol='circle',
line=dict(
color='rgb(104, 74, 114)',
width=1
),
opacity=0.9
)
)
data = [trace1, trace2, trace3]
## Layout settings
layout = go.Layout(
scene = dict(
xaxis = dict(
title= 'sepal length (cm)'),
yaxis = dict(
title= 'petal length (cm)'),
zaxis = dict(
title= 'petal width (cm)'),),
)
fig = go.Figure(data=data, layout=layout)
plot(fig)
这将为您提供以下交互式图表，如下图所示：
Interactive plotting of petal length, petal width, and sepal length
使用这些图表，您可以更好地理解数据并为建模做准备。
超参数
超参数可以被视为确定模型各种属性之一的高级参数，例如复杂性，训练行为和学习率。 这些参数自然与模型参数有所不同，因为它们需要在训练开始之前设置。
例如，K 均值或 K 最近邻中的 k 是这些算法的超参数。 K 均值中的 K 表示要找到的聚类数，K 最近邻中的 k 表示用于进行预测的最近记录数。
调整超参数是任何机器学习项目中提高预测性能的关键步骤。 有多种调优技术，例如网格搜索，随机搜索和贝叶斯优化，但这些技术不在本章范围之内。
让我们通过以下屏幕快照快速浏览 scikit-learn 库中的 K 均值算法参数：
如您所见，有许多参数可以使用，并且至少应查看算法的函数签名，以查看运行算法之前的选项。
让我们一起玩吧。 基线模型将与样本数据一起使用，几乎具有默认设置，如下所示：
from sklearn.datasets.samples_generator import make_blobs
X, y = make_blobs(n_samples=20, centers=3, n_features=3, random_state=42)
k_means = KMeans(n_clusters=3)
y_hat = k_means.fit_predict(X)
y_hat保留集群的成员身份信息，这与原始标签相同，如您在此处看到的：
y_hat
## array([0, 2, 1, 1, 1, 0, 2, 0, 0, 0, 2, 0, 1, 2, 1, 2, 2, 1, 0, 1],
dtype=int32)
y
## array([0, 2, 1, 1, 1, 0, 2, 0, 0, 0, 2, 0, 1, 2, 1, 2, 2, 1, 0, 1])
您可以使用不同的选项来查看它如何影响训练和预测。
损失函数
损失函数通过测量误差来帮助算法在训练过程中更新模型参数，这是预测性能的指标。 损失函数通常表示为：
其中L测量预测值与实际值之间的差。 在训练过程中，此误差被最小化。 不同的算法具有不同的损失函数，迭代次数将取决于收敛条件。
例如，K 均值的损失函数使点与最接近的簇均值之间的平方距离最小化，如下所示：
您将在以下部分中看到详细的实现。
为单个变量实现我们的算法
让我们为单个变量实现 K 均值算法。 您将从一维向量开始，该向量具有 20 条记录，如下所示：
data = [1,2,3,2,1,3,9,8,11,12,10,11,14,25,26,24,30,22,24,27]
trace1 = go.Scatter(
x=data,
y=[0 for x in data],
mode='markers',
name='Data',
marker=dict(
size=12
)
)
layout = go.Layout(
title='1D vector',
)
traces = [trace1]
fig = go.Figure(data=traces, layout=layout)
plot(fig)
如下图所示，将输出以下图表：
我们的目标是找到在数据中可见的3簇。 为了开始实施 K 均值算法，您需要通过选择随机索引来初始化集群中心，如下所示：
n_clusters = 3
c_centers = np.random.choice(X, n_clusters)
print(c_centers)
## [ 1 22 26]
接下来，您需要计算每个点与聚类中心之间的距离，因此请使用以下代码：
deltas = np.array([np.abs(point - c_centers) for point in X])
deltas
array([[ 7, 26, 10],
[ 6, 25, 9],
[ 5, 24, 8],
[ 6, 25, 9],
[ 7, 26, 10],
[ 5, 24, 8],
[ 1, 18, 2],
[ 0, 19, 3],
[ 3, 16, 0],
[ 4, 15, 1],
[ 2, 17, 1],
[ 3, 16, 0],
[ 6, 13, 3],
[17, 2, 14],
[18, 1, 15],
[16, 3, 13],
[22, 3, 19],
[14, 5, 11],
[16, 3, 13],
[19, 0, 16]])
现在，可以使用以下代码来计算集群成员资格：
deltas.argmin(1)
## array([0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1])
现在，您需要使用以下代码来计算记录和群集中心之间的距离：
c_centers = np.array([X[np.where(deltas.argmin(1) == i)[0]].mean() for i in range(3)])
print(c_centers)
## [ 3.625 25.42857143 11.6 ]
这是一次迭代； 您可以继续计算新的群集中心，直到没有任何改善为止。
您可以编写一个函数来包装所有这些功能，如下所示：
def Kmeans_1D(X, n_clusters, random_seed=442):
# Randomly choose random indexes as cluster centers
rng = np.random.RandomState(random_seed)
i = rng.permutation(X.shape[0])[:n_clusters]
c_centers = X[i]
# Calculate distances between each point and cluster centers
deltas = np.array([np.abs(point - c_centers) for point in X])
# Get labels for each point
labels = deltas.argmin(1)
while True:
# Calculate mean of each cluster
new_c_centers = np.array([X[np.where(deltas.argmin(1) == i)[0]].mean() for i in range(n_clusters)])
# Calculate distances again
deltas = np.array([np.abs(point - new_c_centers) for point in X])
# Get new labels for each point
labels = deltas.argmin(1)
# If there's no change in centers, exit
if np.all(c_centers == new_c_centers):
break
c_centers = new_c_centers
return c_centers, labels
c_centers, labels = Kmeans_1D(X, 3)
print(c_centers, labels)
## [11.16666667 25.42857143 2.85714286] [2 2 2 2 2 2 0 0 0 0 0 0 0 1 1 1 1 1 1 1]
让我们使用以下代码绘制集群中心的图表：
trace1 = go.Scatter(
x=X,
y=[0 for num in X],
mode='markers',
name='Data',
marker=dict(
size=12
)
)
trace2 = go.Scatter(
x = c_centers,
y = [0 for num in X],
mode='markers',
name = 'Cluster centers',
marker = dict(
size=12,
color = ('rgb(122, 296, 167)')
)
)
layout = go.Layout(
title='1D vector',
)
traces = [trace1, trace2]
fig = go.Figure(data=traces, layout=layout)
plot(fig)
看下图。 给定的代码将输出以下内容：
您可以清楚地看到，可以为每个元素分配一个聚类中心。
修改我们的算法
现在，您已经了解了单个变量的 K 均值内部，可以将该实现扩展到多个变量，并将其应用于更实际的数据集。
本节中使用的数据集来自 UCI 机器学习存储库，它包括批发分销商的客户信息。 有 440 个具有八项特征的客户。 在下面的列表中，前六个特征与相应产品的年度支出相关，第七个特征显示了购买该产品的渠道，第八个特征显示了该地区：
FRESHMILKGROCERYFROZENDETERGENTS_PAPERDELICATESSENCHANNELREGION
首先下载数据集并将其读取为numpy数组：
from numpy import genfromtxt
wholesales_data = genfromtxt('Wholesale customers data.csv', delimiter=',', skip_header=1)
您可以快速查看数据。 这里是：
print(wholesales_data[:5])
[[2.0000e+00 3.0000e+00 1.2669e+04 9.6560e+03 7.5610e+03 2.1400e+02
2.6740e+03 1.3380e+03]
[2.0000e+00 3.0000e+00 7.0570e+03 9.8100e+03 9.5680e+03 1.7620e+03
3.2930e+03 1.7760e+03]
[2.0000e+00 3.0000e+00 6.3530e+03 8.8080e+03 7.6840e+03 2.4050e+03
3.5160e+03 7.8440e+03]
[1.0000e+00 3.0000e+00 1.3265e+04 1.1960e+03 4.2210e+03 6.4040e+03
5.0700e+02 1.7880e+03]
[2.0000e+00 3.0000e+00 2.2615e+04 5.4100e+03 7.1980e+03 3.9150e+03
1.7770e+03 5.1850e+03]]
选中shape将显示行数和变量数，如下所示：
wholesales_data.shape
## (440, 8)
该数据集具有440个记录，每个记录都有8个特征。
通过使用以下代码来规范化数据集是一个好主意：
wholesales_data_norm = wholesales_data / np.linalg.norm(wholesales_data)
print(wholesales_data_norm[:5])
[[ 1\. 0\. 0.30168043 1.06571214 0.32995207 -0.46657183
0.50678671 0.2638102 ]
[ 1\. 0\. -0.1048095 1.09293385 0.56599336 0.08392603
0.67567015 0.5740085 ]
[ 1\. 0\. -0.15580183 0.91581599 0.34441798 0.3125889
0.73651183 4.87145892]
[ 0\. 0\. 0.34485007 -0.42971408 -0.06286202 1.73470839
-0.08444172 0.58250708]
[ 1\. 0\. 1.02209184 0.3151708 0.28726 0.84957326
0.26205579 2.98831445]]
您可以使用以下代码将数据集读取到pandas.DataFrame中：
import pandas as pd
df = pd.DataFrame(wholesales_data_norm,
columns=['Channel',
'Region',
'Fresh',
'Milk',
'Grocery',
'Frozen',
'Detergents_Paper',
'Delicatessen'])
df.head(10)
让我们创建一个scatterplot矩阵以更仔细地查看数据集。 看一下代码：
fig = ff.create_scatterplotmatrix(df, diag='histogram', size=7, height=1200, width=1200)
plot(fig)
您还可以通过运行以下命令来检查要素之间的相关性：
df.corr()
这将为您提供一个关联表，如下所示：
您还可以使用seaborn热图，如下所示：
import seaborn as sns; sns.set()
ax = sns.heatmap(df.corr(), annot=True)
Correlations between the features
您可以看到特征之间存在一些很强的相关性，例如Grocery和Detergents_Paper之间的相关性。
让我们使用以下代码绘制Grocery，Detergents_Paper和Milk三个特征：
## Creating data for the plotly
trace1 = go.Scatter3d(
# Extracting data based on label
x=df['Grocery'],
y=df['Detergents_Paper'],
z=df['Milk'],
mode='markers',
marker=dict(
size=12,
line=dict(
color='rgba(217, 217, 217, 0.14)',
width=0.5
),
opacity=0.8
)
)
## Layout settings
layout = go.Layout(
scene = dict(
xaxis = dict(
title= 'Grocery'),
yaxis = dict(
title= 'Detergents_Paper'),
zaxis = dict(
title= 'Milk'),),
)
data = [trace1]
fig = dict(data=data, layout=layout)
plot(fig)
现在，您将继续扩展已为更高维度实现的 K 均值算法。 首先，您可以从数据集中删除Channel和Region，如下所示：
df = df[[col for col in df.columns if col not in ['Channel', 'Region']]]
df.head(10)
在实现方面，您还可以使用np.linalg.norm来计算距离，这实际上取决于您使用哪种距离函数。 另一个替代方法是scipy.spatial中的distance.euclidean，如下所示：
def Kmeans_nD(X, n_clusters, random_seed=442):
# Randomly choose random indexes as cluster centers
rng = np.random.RandomState(random_seed)
i = rng.permutation(X.shape[0])[:n_clusters]
c_centers = X[i]
# Calculate distances between each point and cluster centers
deltas = np.array([[np.linalg.norm(i - c) for c in c_centers] for i in X])
# Get labels for each point
labels = deltas.argmin(1)
while True:
# Calculate mean of each cluster
new_c_centers = np.array([X[np.where(deltas.argmin(1) == i)[0]].mean(axis=0) for i in range(n_clusters)])
# Calculate distances again
deltas = np.array([[np.linalg.norm(i - c) for c in new_c_centers] for i in X])
# Get new labels for each point
labels = deltas.argmin(1)
# If there's no change in centers, exit
if np.array_equal(c_centers, new_c_centers):
break
c_centers = new_c_centers
return c_centers, labels
Grocery和Detergents_Paper将用于聚类，并且k将设置为3。 通常，应使用外观检查或弯头方法确定k，如下所示：
centers, labels = Kmeans_nD(df[['Grocery', 'Detergents_Paper']].values, 3)
现在，您可以使用以下命令在数据集中再添加一列：
df['labels'] = labels
您可以使用以下代码来首先可视化结果，以查看结果是否有意义：
## Creating data for the plotly
trace1 = go.Scatter(
# Extracting data based on label
x=df[df['labels'] == 0]['Grocery'],
y=df[df['labels'] == 0]['Detergents_Paper'],
mode='markers',
name='clust_1',
marker=dict(
size=12,
line=dict(
color='rgba(217, 217, 217, 0.14)',
width=0.5
),
opacity=0.8
)
)
trace2 = go.Scatter(
# Extracting data based on label
x=df[df['labels'] == 1]['Grocery'],
y=df[df['labels'] == 1]['Detergents_Paper'],
mode='markers',
name='clust_2',
marker=dict(
color='rgb(#3742fa)',
size=12,
symbol='circle',
line=dict(
color='rgb(204, 204, 204)',
width=1
),
opacity=0.9
)
)
trace3 = go.Scatter(
# Extracting data based on label
x=df[df['labels'] == 2]['Grocery'],
y=df[df['labels'] == 2]['Detergents_Paper'],
mode='markers',
name='clust_3',
marker=dict(
color='rgb(#ff4757)',
size=12,
symbol='circle',
line=dict(
color='rgb(104, 74, 114)',
width=1
),
opacity=0.9
)
)
data = [trace1, trace2, trace3]
## Layout settings
layout = go.Layout(
scene = dict(
xaxis = dict(
title= 'Grocery'),
yaxis = dict(title= 'Detergents_Paper'),
)
)
fig = go.Figure(data=data, layout=layout)
plot(fig)
Plotting the clusters
乍一看，集群看起来很合理，并将最终取决于领域知识支持的解释。
您可以使用以下代码轻松查看每个集群每个特征的平均支出：
df.groupby('labels').mean()
从这个简单的集群中可以看出，第三集群在Milk，Grocery和Detergents_Paper上的支出最高。 第二个集群的支出较低，第一个集群倾向于Milk，Grocery和Detergents_Paper，因此k=2也可以使用。
总结
在本章中，您学习了无监督学习的基础知识，并使用 K 均值算法进行聚类。
有许多聚类算法显示不同的行为。 当涉及到无监督学习算法时，可视化是关键，并且您已经看到了几种不同的方式来可视化和检查数据集。
在下一章中，您将学习 NumPy 常用的其他库，例如 SciPy，Pandas 和 scikit-learn。 这些都是从业者工具包中的重要库，它们相互补充。 您会发现自己将这些库与 NumPy 一起使用，因为每个库都会使某些任务变得更容易。 因此，重要的是要了解有关 Python 数据科学堆栈的更多信息。
六、NumPy，SciPy，Pandas 和 Scikit-Learn
到目前为止，您应该能够使用 NumPy 编写小型实现。 在整个章节中，我们旨在提供使用其他库的示例，在本章中，我们应退后一步，看看可以与 NumPy 一起用于项目的周围库。
本章将介绍其他 Python 库如何对 NumPy 进行补充。 我们将研究以下主题：
- NumPy 和 SciPy
- NumPy 和 Pandas
- SciPy 和 Scikit-learn
NumPy 和 SciPy
到目前为止，您已经看到了许多有关 NumPy 用法的示例，而只有少数 SciPy。 NumPy 具有数组数据类型，它允许您执行各种数组操作，例如排序和整形。
NumPy 具有一些数值算法，可用于执行诸如计算范数，特征值和特征向量之类的任务。 但是，如果数值算法是您的重点，则理想情况下应使用 SciPy，因为它包含更全面的算法集以及最新版本的算法。 SciPy 有许多有用的子程序包可用于某些类型的分析。
以下列表将使您对子软件包有一个整体的了解：
Cluster：此子程序包包含聚类算法。 它具有两个子模块vq和hierarchy。vq模块提供用于 K 均值聚类的函数。 层次结构模块包括用于层次结构聚类的函数。Fftpack：此子程序包包含用于快速傅立叶变换的函数和算法，以及差分和伪差分算符。Interpolate：此子程序包提供用于单变量和多变量插值的函数：1D 和 2D 样条曲线。Linalg：此子程序包提供用于线性代数的函数和算法，例如matrix运算和函数，特征值和-向量计算，矩阵分解，矩阵方程求解器和特殊矩阵。Ndimage：此子程序包提供用于多维图像处理的函数和算法，例如滤镜，插值，测量和形态。Optimize：此子程序包提供函数和算法，用于函数局部和全局优化，函数拟合，求根和线性编程。Signal：此子程序包提供信号处理的函数和算法，例如卷积，B 样条，滤波，连续和离散时间线性系统，波形，小波和频谱分析。Stats：此子程序包提供概率分布，例如连续分布，多元分布和离散分布，以及可以找到均值，众数，方差，偏度，峰度和相关系数的统计函数。
让我们来看一下其中的一个子包。 以下代码显示了用于群集分析的cluster程序包：
Scipy.cluster
%matplotlib inline
import matplotlib.pyplot as plt
## Import ndimage to read the image
from scipy import ndimage
## Import cluster for clustering algorithms
from scipy import cluster
## Read the image
image = ndimage.imread("cluster_test_image.jpg")
## Image is 1000x1000 pixels and it has 3 channels.
image.shape
(1000, 1000, 3)
这将为您提供以下输出：
array([[[30, 30, 30],
[16, 16, 16],
[14, 14, 14],
...,
[14, 14, 14],
[16, 16, 16],
[29, 29, 29]],
[[13, 13, 13],
[ 0, 0, 0],
[ 0, 0, 0],
...,
[ 0, 0, 0],
[ 0, 0, 0],
[12, 12, 12]],
[[16, 16, 16],
[ 3, 3, 3],
[ 1, 1, 1],
...,
[ 0, 0, 0],
[ 2, 2, 2],
[16, 16, 16]],
...,
[[17, 17, 17],
[ 3, 3, 3],
[ 1, 1, 1],
...,
[34, 26, 39],
[27, 21, 33],
[59, 55, 69]],
[[15, 15, 15],
[ 2, 2, 2],
[ 0, 0, 0],
...,
[37, 31, 43],
[34, 28, 42],
[60, 56, 71]],
[[33, 33, 33],
[20, 20, 20],
[17, 17, 17],
...,
[55, 49, 63],
[47, 43, 57],
[65, 61, 76]]], dtype=uint8)
在这里，您可以看到该图：
plt.figure(figsize = (15,8))
plt.imshow(image)
您可以从前面的代码块中获得以下图表：
使用以下代码将图像数组转换为二维数据集：
x, y, z = image.shape
image_2d = image.reshape(x*y, z).astype(float)
image_2d.shape
(1000000, 3)
image_2d
array([[30., 30., 30.],
[16., 16., 16.],
[14., 14., 14.],
...,
[55., 49., 63.],
[47., 43., 57.],
[65., 61., 76.]])
## kmeans will return cluster centers and the distortion
cluster_centers, distortion = cluster.vq.kmeans(image_2d, k_or_guess=2)
print(cluster_centers, distortion)
[[179.28653454 179.30176248 179.44142117]
[ 3.75308484 3.83491111 4.49236356]] 26.87835069294931
image_2d_labeled = image_2d.copy()
labels = []
from scipy.spatial.distance import euclidean
import numpy as np
for i in range(image_2d.shape[0]):
distances = [euclidean(image_2d[i], center) for center in cluster_centers]
labels.append(np.argmin(distances))
plt.figure(figsize = (15,8))
plt.imshow(cluster_centers[labels].reshape(x, y, z))
您从前面的代码中获得以下输出：
SciPy 和 NumPy 和线性回归
您已经了解了如何使用 NumPy 从头开始编写线性回归算法。Scipy.stats模块具有linregress函数，用于计算斜率，截距，相关系数（r 值），两侧 p 值以及标准差估计，如下所示：
from sklearn import datasets
%matplotlib inline
import matplotlib.pyplot as plt
## Boston House Prices dataset
boston = datasets.load_boston()
x = boston.data
y = boston.target
boston.feature_names
array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
'TAX', 'PTRATIO', 'B', 'LSTAT'], dtype='<U7')
x.shape
(506, 13)
y.shape
(506,)
## We will consider "lower status of population" as independent variable for its importance
lstat = x[0:,-1]
lstat.shape
(506,)
from scipy import stats
slope, intercept, r_value, p_value, std_err = stats.linregress(lstat, y)
print(slope, intercept, r_value, p_value, std_err)
-0.9500493537579909 34.55384087938311 -0.737662726174015 5.081103394387796e-88 0.03873341621263942
print("r-squared:", r_value**2)
r-squared: 0.5441462975864798
plt.plot(lstat, y, 'o', label='original data')
plt.plot(lstat, intercept + slope*lstat, 'r', label='fitted line')
plt.legend()
plt.show()
我们从前面代码的输出中获得以下图表，如下图所示：
您还可以查看平均房间数与房价之间的关系。 以下代码块打印出性能指标：
rm = x[0:,5]
slope, intercept, r_value, p_value, std_err = stats.linregress(rm, y)
print(slope, intercept, r_value, p_value, std_err)
print("r-squared:", r_value**2)
## 9.102108981180308 -34.670620776438554 0.6953599470715394 2.48722887100781e-74 0.4190265601213402
## r-squared: 0.483525455991334
以下代码块绘制了拟合线：
plt.plot(rm, y, 'o', label='original data')
plt.plot(rm, intercept + slope*rm, 'r', label='fitted line')
plt.legend()
plt.show()
我们从前面的代码中获得以下输出，如下图所示：
NumPy 和 Pandas
考虑一下时，NumPy 是一个相当低级的数组操作库，大多数其他 Python 库都写在它的顶部。
这些库之一是pandas，它是一个高级数据处理库。 浏览数据集时，通常会执行诸如计算描述性统计数据，按特定特征分组以及合并之类的操作。pandas库具有许多友好的函数来执行这些各种有用的操作。
在此示例中，我们使用糖尿病数据集。sklearn.datasets中的糖尿病数据集使用零均值和 L2 范数标准化。
该数据集包含 442 条记录，这些记录具有 10 个特征：年龄，性别，体重指数，平均血压和 6 个血清测量值。
目标代表采取这些基准措施后的疾病进展。 您可以在 web 和相关论文 中查看数据描述。
我们从操作开始，如下所示：
import pandas as pd
from sklearn import datasets
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
diabetes = datasets.load_diabetes()
df = pd.DataFrame(diabetes.data, columns=diabetes.feature_names)
diabetes.feature_names
['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6']
df.head(10)
我们从前面的代码中获得以下输出，如下表所示：
这段代码向您展示了如何在数据帧中创建目标列：
df['Target'] = diabetes.target
df.head(10)
Pandas 帮助我们轻松地处理表格数据，并通过各种辅助方法和可视化支持我们的分析。 看一下代码：
## Descriptive statistics
df.describe()
我们从前面的代码中获得以下输出，如下表所示：
通过使用以下代码行，看看目标的分布方式：
plt.hist(df['Target'])
下图显示了上一行的输出：
您可以看到目标变量向右倾斜。 看一下这段代码：
## Since 'sex' is categorical, excluding it from numerical columns
numeric_cols = [col for col in df.columns if col != 'sex']
numeric_cols
## ['age', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6', 'Target']
## You can have a look at variable distributions individually, but there's a better way
df[numeric_cols].hist(figsize=(20, 20), bins=30, xlabelsize=12, ylabelsize=12)
## You can also choose create dataframes for numerical and categorical variables
前一个代码块的输出：
Feature distributions
您可以检查其中一些特征的分布，并确定其中哪些看起来相似。 对于此示例，特征s1，s2和s6似乎具有相似的分布，如从此代码中可以看到的：
## corr method will give you the correlation between features
df[numeric_cols].corr()
下图显示了上一行的输出：
您可以使用heatmap更好地表示这种关系，如下所示：
plt.figure(figsize=(15, 15))
sns.heatmap(df[numeric_cols].corr(), annot=True)
下图是由前面的代码块生成的热图：
Correlations heatmap
您还可以通过以下方式过滤相关性：
plt.figure(figsize=(18, 15))
sns.heatmap(df[numeric_cols].corr()
[(df[numeric_cols].corr() >= 0.3) & (df[numeric_cols].corr() <= 0.5)],
annot=True)
此图显示了过滤后的相关性：
Filtered correlations heatmap
您还可以使用其他有用的可视化来检查统计关系，如下所示：
fig, ax = plt.subplots(3, 3, figsize = (18, 12))
for i, ax in enumerate(fig.axes):
if i < 9:
sns.regplot(x=df[numeric_cols[i]],y='Target', data=df, ax=ax)
该图显示了前面代码的以下输出：
Regression Plots
您可以看到，使用pandas使探索性数据分析相对简单。 使用pandas，您可以检查特征及其关系。
Pandas 和股票价格的定量建模
pandas最初是为在金融数据集中使用而编写的，它包含许多用于处理时间序列数据的便捷函数。 在本节中，您将看到如何使用pandas库处理股票价格序列。
您将使用quandl Python 库获取公司的财务数据。 看一下这段代码：
import quandl msft = quandl.get('WIKI/MSFT')
msft.columns
## Index(['Open', 'High', 'Low', 'Close', 'Volume', 'Ex-Dividend', 'Split Ratio', 'Adj. Open', 'Adj. High', 'Adj. Low', 'Adj. Close', 'Adj. Volume'], dtype='object')
msft.tail()
下表显示了msft.tail()的输出：
让我们通过导入以下设置来自定义绘图：
matplotlib.font_manager as font_manager font_path = '/Library/Fonts/Cochin.ttc' font_prop = font_manager.FontProperties(fname=font_path, size=24) axis_font = {'fontname':'Arial', 'size':'18'} title_font = {'fontname':'Arial', 'size':'22', 'color':'black', 'weight':'normal', 'verticalalignment':'bottom'} plt.figure(figsize=(10, 8)) plt.plot(msft['Adj. Close'], label='Adj. Close') plt.xticks(fontsize=22) plt.yticks(fontsize=22) plt.xlabel("Date", **axis_font) plt.ylabel("Adj. Close", **axis_font) plt.title("MSFT", **title_font) plt.legend(loc='upper left', prop=font_prop, numpoints=1) plt.show()
该图显示了来自先前设置的图：
您可以使用以下代码来计算每日变化：
msft['Daily Pct. Change'] = (msft['Adj. Close'] - msft['Adj. Open']) / msft['Adj. Open']
msft.tail(10)
下图显示了msft.tail(10)的输出：
尝试使用此每日收益的直方图：
plt.figure(figsize=(22, 8))
plt.hist(msft['Daily Pct. Change'], bins=100)
如下图所示，它将为您提供以下图表：
Distribution of daily returns
返回的分布具有较长的尾巴，尤其是在负侧，这是财务分析中的已知现象。 它产生的风险称为尾部风险，它与市场收益服从正态分布的假设相矛盾。 这基本上告诉您，极端事件发生的可能性比更正态分布的可能性更大。
在可视化方面，使它们具有交互性很有帮助。 为此，plotly提供了一个很好的替代当前绘图库的方法，如下所示：
import plotly.plotly as py
import plotly.graph_objs as go
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot init_notebook_mode(connected=True)
from datetime import datetime
import pandas_datareader.data as web
import quandl
msft = quandl.get('WIKI/MSFT')
msft['Daily Pct. Change'] = (msft['Adj. Close'] - msft['Adj. Open']) / msft['Adj. Open']
data = [go.Scatter(x=msft.index, y=msft['Adj. Close'])]
plot(data)
我们从前面的代码中获得以下图表，如下图所示：
您可以创建开盘-最高-最低-收盘价（OHLC）的图表，其中每个日期都有 4 个不同的价格值，它们是开，高，低和关。 看一下这段代码：
charts trace = go.Ohlc(x=msft.index, open=msft['Adj. Open'], high=msft['Adj. High'], low=msft['Adj. Low'], close=msft['Adj. Close'])
data = [trace]
plot(data)
该图显示了先前代码的图：
您可以通过在图表上选择自定义范围来检查特定区域。 看一下这个图：
同样，您可以使用以下代码创建Candlestick图表：
trace = go.Candlestick(x=msft.index, open=msft['Adj. Open'], high=msft['Adj. High'], low=msft['Adj. Low'], close=msft['Adj. Close'])
data = [trace]
plot(data)
下图显示了此代码的输出：
您也可以在Candlestick图表中选择特定范围。 看一下这个图：
分布图如下：
import plotly.figure_factory as ff
fig = ff.create_distplot([msft['Daily Pct. Change'].values], ['MSFT Daily Returns'], show_hist=False)
plot(fig)
下图显示了前面代码的输出：
我们可以创建三个移动平均线，如下所示：
msft['200MA'] = msft['Adj. Close'].rolling(window=200).mean()
msft['100MA'] = msft['Adj. Close'].rolling(window=100).mean()
msft['50MA'] = msft['Adj. Close'].rolling(window=50).mean()
msft.tail(10)
下表显示了msft.tail(10)的输出：
根据切片的数据，将包括最近 2,000 天。 看一下这段代码：
trace_adjclose = go.Scatter( x=msft[-2000:].index, y=msft[-2000:]['Adj. Close'], name = "Adj. Close", line = dict(color = '#000000'), opacity = 0.8)
trace_200 = go.Scatter( x=msft[-2000:].index, y=msft[-2000:]['200MA'], name = "200MA", line = dict(color = '#FF0000'), opacity = 0.8)
trace_100 = go.Scatter( x=msft[-2000:].index, y=msft[-2000:]['100MA'], name = "100MA", line = dict(color = '#0000FF'), opacity = 0.8)
trace_50 = go.Scatter( x=msft[-2000:].index, y=msft[-2000:]['50MA'], name = "50MA", line = dict(color = '#FF00FF'), opacity = 0.8)
data = [trace_adjclose, trace_200, trace_100, trace_50]
layout = dict( title = "MSFT Moving Averages: 200, 100, 50 days", )
fig = dict(data=data, layout=layout)
plot(fig)
下图显示了前面代码块中的图：
移动平均线用于监视金融市场的趋势。 在此示例中，有三个移动平均线，每个移动线均具有不同的周期。 您可以设置分析天数，以进行短期，中期和长期趋势监视。
当您开始使用财务时间序列时，您会很快意识到您需要基于不同期间的汇总，并且在pandas中创建这些汇总非常容易。 以下代码段将通过计算平均值每月汇总记录：
msft_monthly = msft.resample('M').mean()
msft_monthly.tail(10)
下图显示了msft_monthly.tail(10)的输出：
这是一个简单的时间序列图：
data = [go.Scatter(x=msft_monthly[-24:].index, y = msft_monthly[-24:]['Adj. Close'])]
plot(data)
如下图所示，这将为您提供以下图表：
在检查要素之间的关系时，可以使用我们在前面的示例中已经看到的相关矩阵。 在时间序列中，从业者对自相关感兴趣，自相关显示了时间序列与其自身的相关性。 例如，理想情况下，您希望时间序列中出现周期性的峰值，以显示您的季节性。 通过使用以下代码，让我们看看每日百分比变化是否有任何明显的峰值：
plt.figure(figsize=(22, 14)) pd.plotting.autocorrelation_plot(msft_monthly['Daily Pct. Change'])
我们从前面的代码中获得以下图表，如下图所示：
Monthly autocorrelation plot
在这个系列中没有有意义的显着滞后，但是如果您使用宏观经济变量（例如 GDP，通货膨胀率和失业水平）进行尝试，则可能会看到显着的季度或年度峰值。
SciPy 和 Scikit-learn
Scikit-learn 是用于机器学习的 SciKit 库之一，它建立在 SciPy 之上。 您可以使用它执行回归分析，就像在前几章中使用 scikit-learn 库所做的那样。 看一下这段代码：
from sklearn import datasets, linear_model
from sklearn.metrics import mean_squared_error, r2_score
diabetes = datasets.load_diabetes()
linreg = linear_model.LinearRegression()
linreg.fit(diabetes.data, diabetes.target)
## You can inspect the results by looking at evaluation metrics
print('Coeff.: n', linreg.coef_)
print("MSE: {}".format(mean_squared_error(diabetes.target, linreg.predict(diabetes.data)))) print('Variance Score: {}'.format(r2_score(diabetes.target, linreg.predict(diabetes.data))))
scikit-learn 和住房数据的 K 均值聚类
在本节中，我们将使用 scikit-learn 的 K 均值算法对房屋数据进行聚类，如下所示：
from sklearn.cluster import KMeans from sklearn.datasets import load_boston boston = load_boston()
## As previously, you have implemented the KMeans from scracth and in this example, you use sklearns API k_means = KMeans(n_clusters=3) # Training k_means.fit(boston.data)
KMeans(algorithm='auto', copy_x=True, init='K 均值++', max_iter=300, n_clusters=3, n_init=10, n_jobs=1, precompute_distances='auto', random_state=None, tol=0.0001, verbose=0)
print(k_means.labels_)
前面代码的输出如下：
[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 0 2 2
2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 0 0 0 0 0 0 0 0 0 0 0 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 2
2 0 2 2 2 2 0 2 2 2 0 0 0 0 2 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2
2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1]
您可以使用以下代码行检查群集中心：
print(k_means.cluster_centers_)
这是控制台的输出：
[[ 1.49558803e+01 -5.32907052e-15 1.79268421e+01 2.63157895e-02
6.73710526e-01 6.06550000e+00 8.99052632e+01 1.99442895e+00
2.25000000e+01 6.44736842e+02 1.99289474e+01 5.77863158e+01
2.04486842e+01]
[ 3.74992678e-01 1.57103825e+01 8.35953552e+00 7.10382514e-02
5.09862568e-01 6.39165301e+00 6.04133880e+01 4.46074481e+00
4.45081967e+00 3.11232240e+02 1.78177596e+01 3.83489809e+02
1.03886612e+01]
[ 1.09105113e+01 5.32907052e-15 1.85725490e+01 7.84313725e-02
6.71225490e-01 5.98226471e+00 8.99137255e+01 2.07716373e+00
2.30196078e+01 6.68205882e+02 2.01950980e+01 3.71803039e+02
1.78740196e+01]]
关于聚类算法的评估，通常将使用诸如轮廓分析或弯头方法之类的技术来评估聚类的质量并确定正确的超参数（例如 K 均值）。 使用 scikit-learn 提供的简单 API，您还将发现这种分析易于执行。 强烈建议您通过实践从这些示例中学到的知识为基础，这将提高您的知识和技能。
总结
在本章中，您使用各种示例（主要用于机器学习任务）练习了 NumPy，SciPy，Pandas 和 scikit-learn。 使用 Python 数据科学库时，通常有不止一种执行给定任务的方法，而且通常有助于了解不止一种方法。
您可以使用替代方法以获得更好的实现，也可以出于比较的目的。 在为给定任务尝试不同的方法时，您可能会找到不同的选项，这些选项将允许您进一步自定义实现，或者只是观察到一些性能改进。
本章的目的是向您展示这些不同的选项，以及 Python 语言由于其丰富的分析库生态系统而具有的灵活性。 在下一章中，您将了解有关 NumPy 内部的更多信息，例如 numpy 如何管理数据结构和内存，代码概要分析以及有效编程的技巧。
七、高级 NumPy
许多库都具有易于使用的 API。 您需要做的就是调用提供的 API 函数，库将为您处理其余的函数。 幕后发生的事情与您无关，您只对输出感兴趣。 在大多数情况下，这很好，但是，至少要了解所使用库的基本内部结构很重要。 了解内部结构将帮助您掌握代码的最新动态以及开发应用时应避免的危险信号。
在本章中，将回顾 NumPy 的内部结构，例如 NumPy 的类型层次结构和内存使用情况。 在本章的最后，您将了解用于逐行检查程序的代码配置文件。
NumPy 内部
如您在前几章中所见，NumPy 数组使数值计算变得高效，其 API 直观且易于使用。 NumPy 数组也是其他科学图书馆的核心，因为其中许多库都是基于 NumPy 数组构建的。
为了编写更好，更高效的代码，您需要了解数据处理的内部。 NumPy 数组及其元数据位于数据缓冲区中，该缓冲区是带有某些数据项的专用内存块。
NumPy 如何管理内存？
初始化 NumPy 数组后，其元数据和数据将存储在随机存取存储器（RAM）中分配的存储位置上。
import numpy as np
array_x = np.array([100.12, 120.23, 130.91])
首先，Python 是一种动态类型化的语言； 不需要显式声明变量类型，例如int或double。 可以推断变量类型，在这种情况下，您希望array_x的数据类型为np.float64：
print(array_x.dtype)
float64
使用numpy库而不是 Python 的优势在于numpy支持许多不同的数值数据类型，例如bool_，int_，intc，intp，int8，int16，int32，int64，uint8，uint16，uint32，uint64，float_，float16，float32，float64，complex_，complex64和complex128。
您可以通过检查sctypes查看这些类型：
np.sctypes
{'complex': [numpy.complex64, numpy.complex128, numpy.complex256],
'float': [numpy.float16, numpy.float32, numpy.float64, numpy.float128],
'int': [numpy.int8, numpy.int16, numpy.int32, numpy.int64],
'others': [bool, object, bytes, str, numpy.void],
'uint': [numpy.uint8, numpy.uint16, numpy.uint32, numpy.uint64]}
下图显示了数据类型树：
您可以通过调用mro方法来检查诸如np.float64之类的数据类型的父类：
np.float64.mro()
[numpy.float64,
numpy.floating,
numpy.inexact,
numpy.number,
numpy.generic,
float,
object]
和np.int64的父类：
np.int64.mro()
[numpy.int64,
numpy.signedinteger,
numpy.integer,
numpy.number,
numpy.generic,
object]
mro方法代表“方法解析顺序”。 为了更好地理解继承，应该首先了解继承的概念。 在可以使用面向对象范例的编程语言中，可以将对象的属性和方法基于之前创建的另一个对象，这就是继承。 在前面的示例中，np.int64保留了np.signedinteger及其之后的属性和行为。
让我们看一个简单的例子：
class First:
def firstmethod(self):
print("Call from First Class, first method.")
class Second:
def secondmethod(self):
print("Call from Second Class, second method.")
class Third(First, Second):
def thirdmethod(self):
print("Call from Third Class, third method.")
这里，有 3 个类别，而First和Second类别是独立的，Third类别是从First和Second继承的。 您可以创建Third类的实例，并使用dir方法检查其内容：
myclass = Third()
dir(myclass)
[...
'__repr__',
'__setattr__',
'__sizeof__',
'__str__',
'__subclasshook__',
'__weakref__',
'firstmethod',
'secondmethod',
'thirdmethod']
dir表示myclass的方法中有firstmethod，secondmethod和thirdmethod。
您可以调用这些方法，如下所示：
myclass.firstmethod()
myclass.secondmethod()
myclass.thirdmethod()
## Call from First Class, first method.
## Call from Second Class, second method.
## Call from Third Class, third method.
现在，让我们向Second类添加firstmethod，看看会发生什么：
class First:
def firstmethod(self):
print("Call from First Class, first method.")
class Second:
def firstmethod(self):
print("Call from Second Class, first method.")
def secondmethod(self):
print("Call from Second Class, second method.")
class Third(First, Second):
def thirdmethod(self):
print("Call from Third Class, third method.")
像以前一样检查方法输出：
myclass = Third()
myclass.firstmethod()
myclass.secondmethod()
myclass.thirdmethod()
## Call from First Class, first method.
## Call from Second Class, second method.
## Call from Third Class, third method.
如您所见，已添加到Second类的方法无效，因为Third类的实例从First类继承了该实例。
您可以按以下方式检查类的mro：
Third.__mro__
这将为您提供以下输出：
(__main__.Third, __main__.First, __main__.Second, object)
这是使用继承机制时解析属性和方法的方式，并且现在您应该或多或少地了解mro的工作方式。 现在，您可以再次查看我们之前拥有的 numpy 数据类型的mro示例。
您可以使用nbytes来查看存储数据类型所需的内存。
首先，让我们看看单个float64的大小：
np.float64(100.12).nbytes
8
np.str_('n').nbytes
4
np.str_('numpy').nbytes
20
array_x具有 3 个float64，其大小将是元素数乘以商品大小，即24，如以下代码段所示：
np.float64(array_x).nbytes
24
例如，如果您需要较低的计算精度，则可以使用np.float32，它将占用float64占用的一半内存：
array_x2 = array_x.astype(np.float32)
array_x2
array([100.12, 120.23, 130.91], dtype=float32)
np.float32(array_x2).nbytes
12
简单来说，8 个字节的内存将保存 1 float64或 2 float32。
Python 的动态性质引入了一种处理数据类型的新方法，因为 Python 应该包含有关其存储的数据的更多信息。 虽然典型的 C 变量将具有有关内存位置的信息，但 Python 变量应具有存储为 C 结构的信息，该结构包含引用计数，对象的类型，对象的大小以及变量本身。
这是提供灵活的环境来处理不同数据类型所必需的。 如果诸如列表之类的数据结构可以容纳不同类型的对象，这是由于该信息对于列表中的每个元素的存储。
但是，由于 NumPy 数组中的数据类型是固定的，由于使用了连续的内存块，因此内存使用效率可能更高。
您可以通过检查 NumPy 数组的__array_interface__属性来查看地址和其他相关信息。
编写此接口是为了允许开发人员共享数组内存和信息：
array_x.__array_interface__
{'data': (140378873611440, False),
'descr': [('', '<f8')],
'shape': (3,),
'strides': None,
'typestr': '<f8',
'version': 3}
__array_interface__是具有 6 个键的 python 字典：
shape的工作方式类似于 NumPy 数组或pandas数据帧的常规shape属性。 它显示每个大小的大小。 由于array_x具有1大小和3元素，因此它是具有3大小的元组。typestr具有3值，第一个显示字节顺序，第二个显示字符代码，其余字符显示字节数。 在此示例中，其值为'<f8'，这表示字节顺序为低位字节序，字符代码为浮点，并且使用的字节数为 8。descr可能会提供有关内存布局的更多详细信息。 默认值为[('', typestr)]。data显示数据的存储位置。 这是一个元组，其中第一个元素显示 NumPy 数组的存储块地址，第二个元素是指示其是否为只读的标志。 在此示例中，内存块地址为140378873611440，它不是只读的。strides指示给定的数组是否为 C 样式的连续内存缓冲区。 在此示例中，None 表示这是 C 样式的连续数组。 否则，它将包含跨步元组，以了解跳转到给定维度中的下一个数组元素所要跳转的位置。 步幅是重要的属性，因为当您使用不同的切片（例如X[::4]）时，步幅将引导数组视图。version表示在此示例中版本号为 3。
以下片段显示了一个简单的示例：
import numpy as np
X = np.array([1,2,3,2,1,3,9,8,11,12,10,11,14,25,26,24,30,22,24,27])
X[::4]
## array([ 1, 1, 11, 14, 30])
这一点很重要，因为当您使用基于现有ndarrays的切片创建新的ndarrays时，可能会降低性能。 让我们看一个简单的例子； 以下代码段创建了 3D ndarray：
nd_1 = np.random.randn(4, 6, 8)
nd_1
## array([[[ 0.64900179, -0.00494884, -0.97565618, -0.78851039],
[ 0.05165607, 0.068948 , 1.54542042, 1.68553396],
[-0.80311258, 0.95298682, -0.85879725, 0.67537715]],
[[ 0.24014811, -0.41894241, -0.00855571, 0.43483418],
[ 0.43001636, -0.75432657, 1.16955535, -0.42471807],
[ 0.6781286 , -1.87876591, 1.02969921, 0.43215107]]])
您可以对其进行切片并创建另一个数组：
nd_2 = nd_1[::, ::2, ::2]
将会选择：
- 首先，第一维的所有项目
- 然后，第二维的每两个项目
- 然后，第三维的每两个项目
它将具有以下数组：
print(nd_2)
[[[ 0.64900179 -0.97565618]
[-0.80311258 -0.85879725]]
[[ 0.24014811 -0.00855571]
[ 0.6781286 1.02969921]]]
您可以看到nd_1和nd_2的内存地址相同：
nd_1.__array_interface__
{'data': (140547049888960, False),
'descr': [('', '<f8')],
'shape': (2, 3, 4),
'strides': None,
'typestr': '<f8',
'version': 3}
nd_2.__array_interface__
{'data': (140547049888960, False),
'descr': [('', '<f8')],
'shape': (2, 2, 2),
'strides': (96, 64, 16),
'typestr': '<f8',
'version': 3}
nd_2大步前进，了解如何沿nd_1数组的不同维度移动。
为了强调这些跨步在数值计算中的作用，下面的示例将为数组维和切片使用更大的大小：
nd_1 = np.random.randn(400, 600)
nd_2 = np.random.randn(400, 600*20)[::, ::20]
nd_1和nd_2具有相同的大小：
print(nd_1.shape, nd_2.shape)
(400, 600) (400, 600)
您可以测量用于计算nd_1和nd_2的数组元素的累积乘积的时间：
%%timeit
np.cumprod(nd_1)
## 802 µs ± 20.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
np.cumprod(nd_2)
## 12 ms ± 71.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
两次操作之间存在明显的时间差； 这是为什么？ 如您所料，nd_2中的步幅过大会导致此问题：
nd_1.__array_interface__
{'data': (4569473024, False),
'descr': [('', '<f8')],
'shape': (400, 600),
'strides': None,
'typestr': '<f8',
'version': 3}
nd_2.__array_interface__
{'data': (4603252736, False),
'descr': [('', '<f8')],
'shape': (400, 600),
'strides': (96000, 160),
'typestr': '<f8',
'version': 3}
从存储器向 CPU 读取数据时，nd_2中存在跨步会导致跳转到不同的存储器位置。 如果将数组元素顺序存储为连续的内存块，那么从时间测量来看，此操作会更快。 步伐越小越好，可以更好地利用 CPU 缓存来提高性能。
有一些变通办法可以缓解与 CPU 缓存相关的问题。 您可以使用的一种库是numexpr库，它是 NumPy 的快速数值表达式求值器。 库使内存使用效率更高，并且还可以使多线程编程受益，以充分利用可用的内核。 您也可以将其与 Intel 的 VML 结合使用以进行进一步的改进。
在下面的示例中，让我们看看它是否对nd_2有帮助：
import numexpr as ne
%%timeit
2 * nd_2 + 48
## 4 ms ± 10.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%%timeit
ne.evaluate("2 * nd_2 + 48")
## 843 µs ± 8.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
您应该使用其他示例进行尝试，以查看性能提升。
如果从头开始索引数组到某个元素，您将看到它具有相同的内存地址：
array_x[:2].__array_interface__['data'][0]
## 140378873611440
array_x[:2].__array_interface__['data'][0] == array_x.__array_interface__['data'][0]
## True
但是，如果您在0以外的位置开始索引，则将为您提供不同的内存地址：
array_x[1:].__array_interface__['data'][0]
## 140378873611448
array_x[1:].__array_interface__['data'][0] == array_x.__array_interface__['data'][0]
## False
ndarray的另一个属性称为标志，它提供有关给定 NumPy 数组的内存布局的信息：
array_f = np.array([[100.12, 120.23, 130.91], [90.45, 110.32, 120.32]])
print(array_f)
## [[100.12 120.23 130.91]
## [ 90.45 110.32 120.32]]
array_f.flags
## C_CONTIGUOUS : True
## F_CONTIGUOUS : False
## OWNDATA : True
## WRITEABLE : True
## ALIGNED : True
## WRITEBACKIFCOPY : False
## UPDATEIFCOPY : False
您可以使用类似于字典的符号或小写的属性名称来获得单个标志：
array_f.flags['C_CONTIGUOUS']
## True
array_f.flags.c_contiguous
## True
让我们看一下每个属性：
C_CONTIGUOUS：C 样式的连续内存的单个块F_CONTIGUOUS：连续内存的单个块，Fortran 风格
您的数据可以使用不同的布局存储在内存中。 这里有 2 种不同的内存布局要考虑：对应于C_CONTIGUOUS的行主要顺序和对应于F_CONTIGUOUS的列主要顺序。
在该示例中，array_f是二维的，array_f的行项目存储在相邻的存储位置中。 类似地，在F_CONTIGUOUS情况下，每列的值存储在相邻的存储位置中。
某些numpy函数将使用参数order将此顺序指示为'C'或'F'。 以下示例显示了具有不同顺序的reshape函数：
np.reshape(array_f, (3, 2), order='C')
## array([[100.12, 120.23],
## [130.91, 90.45],
## [110.32, 120.32]])
np.reshape(array_f, (3, 2), order='F')
## array([[100.12, 110.32],
## [ 90.45, 130.91],
## [120.23, 120.32]])
其余的：
OWNDATA：数组是否与另一个对象共享其内存块或拥有所有权WRITEABLE：False表示它是只读的； 否则可以将该区域写入。ALIGNED：数据是否针对硬件对齐WRITEBACKIFCOPY：该数组是否是另一个数组的副本UPDATEIFCOPY：（不建议使用WRITEBACKIFCOPY）
了解内存管理很重要，因为它会影响性能。 根据您执行计算的方式，计算速度会有所不同。 您可能没有意识到某些计算涉及现有数组的隐式副本，这会减慢计算速度。
以下代码块显示了两个示例，其中第一个不需要复制，而第二个具有隐式复制操作：
shape = (400,400,400)
array_x = np.random.random_sample(shape)
import cProfile
import re
cProfile.run('array_x *= 2')
## 3 function calls in 0.065 seconds
## Ordered by: standard name
## ncalls tottime percall cumtime percall filename:lineno(function)
## 1 0.065 0.065 0.065 0.065 <string>:1(<module>)
## 1 0.000 0.000 0.065 0.065 {built-in method builtins.exec}
## 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}
import cProfile
import re
cProfile.run('array_y = array_x * 2')
## 3 function calls in 0.318 seconds
## Ordered by: standard name
## ncalls tottime percall cumtime percall filename:lineno(function)
## 1 0.318 0.318 0.318 0.318 <string>:1(<module>)
## 1 0.000 0.000 0.318 0.318 {built-in method builtins.exec}
## 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}
首次运行比第二次慢 5 倍。 您需要了解隐式复制操作，并熟悉它在哪种情况下发生。 重塑数组不需要隐式复制，除非它是矩阵转置。
许多数组操作会返回一个新数组以获取结果。 此行为是预期的，但会破坏迭代任务的性能，在迭代任务中，您可能会有数百万或数十亿次迭代。 某些numpy函数具有out参数，该参数创建输出数组，并使用其写入迭代结果。 通过这种方式，您的程序可以更好地管理内存，并且需要更少的时间：
shape_x = (8000,3000)
array_x = np.random.random_sample(shape_x)
%%timeit
np.cumprod(array_x)
## 176 ms ± 2.32 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
output_array的类型和大小应与操作的预期输出相同：
output_array = np.zeros(array_x.shape[0] * array_x.shape[1])
%%timeit
np.cumprod(array_x, out=output_array)
## 86.4 ms ± 1.21 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
对 NumPy 代码进行性能分析以了解性能
有几个有用的库可以监视给定 python 脚本的性能指标。 您已经看到cProfile库的用法。 本节将介绍vprof，它是可视分析器库。
它将为您提供给定 python 程序的运行时统计信息和内存利用率。
第 5 章“使用 NumPy 聚类批发分销商的客户”的一维特征将在此处使用，以下代码段应保存到名为to_be_profiled.py的文件中：
import numpy as np
X = np.array([1,2,3,2,1,3,9,8,11,12,10,11,14,25,26,24,30,22,24,27])
n_clusters = 3
def Kmeans_1D(X, n_clusters, random_seed=442):
# Randomly choose random indexes as cluster centers
rng = np.random.RandomState(random_seed)
i = rng.permutation(X.shape[0])[:n_clusters]
c_centers = X[i]
# Calculate distances between each point and cluster centers
deltas = np.array([np.abs(point - c_centers) for point in X])
# Get labels for each point
labels = deltas.argmin(1)
while True:
# Calculate mean of each cluster
new_c_centers = np.array([X[np.where(deltas.argmin(1) == i)[0]].mean() for i in range(n_clusters)])
# Calculate distances again
deltas = np.array([np.abs(point - new_c_centers) for point in X])
# Get new labels for each point
labels = deltas.argmin(1)
# If there's no change in centers, exit
if np.all(c_centers == new_c_centers):
break
c_centers = new_c_centers
return c_centers, labels
c_centers, labels = Kmeans_1D(X, 3)
print(c_centers, labels)
保存此文件后，您可以使用命令行开始对其进行性能分析。
可以通过 4 种不同的方式配置vprof来获取：
- CPU 火焰图（
vprof -c c to_be_profiled.py） - 内置的探查器统计信息（
vprof -c p to_be_profiled.py） - 运行程序中的行后，Cpython 垃圾收集器跟踪的对象和进程内存的内存图（
vprof -c m to_be_profiled.py） - 用
runtime编码heatmap并为每条执行的行计数（vprof -c h to_be_profiled.py）
这 4 种配置可用于单个源文件或包。 让我们看一下p，m和h配置的输出：
探查器的配置：
$ vprof -c p to_be_profiled.py
Running Profiler...
[10.71428571 25.42857143 2\. ] [2 2 2 2 2 2 0 0 0 0 0 0 0 1 1 1 1 1 1 1]
Starting HTTP server...
一个新的选项卡将在您的浏览器中打开，并显示以下输出：
Time spent for each call
您可以看到文件名，函数名，行号和每次调用所花费的时间。
内存使用情况统计信息的配置：
$ vprof -c m to_be_profiled.py
Running MemoryProfiler...
[10.71428571 25.42857143 2\. ] [2 2 2 2 2 2 0 0 0 0 0 0 0 1 1 1 1 1 1 1]
Starting HTTP server...
新标签页将在您的浏览器中打开，并显示以下输出：
内存使用情况
在左侧，您可以看到内存中的对象，并且图表显示了随着执行的行数增加，内存使用量（以兆字节为单位）。 如果将鼠标悬停在图表上，则每个点将具有以下信息：
- 执行的行
- 行号
- 函数名称
- 文件名称
- 内存使用情况
例如，to_be_profiled.py文件中的第 27 行是计算deltas的下一行：
deltas = np.array([np.abs(point - new_c_centers) for point in X])
它执行了很多次，因为如果您检查图表，这是列表理解。
代码heatmap的配置：
$ vprof -c h to_be_profiled.py
Running CodeHeatmapProfiler...
[10.71428571 25.42857143 2\. ] [2 2 2 2 2 2 0 0 0 0 0 0 0 1 1 1 1 1 1 1]
Starting HTTP server...
“新”标签将在您的浏览器中打开，并显示以下输出：
Heatmap for the lines executed
在左侧，有一个已检查模块的列表，在这种情况下，只有一个文件要检查。
右边是程序中每行的热图。 如果将鼠标悬停在行上，它将为您提供以下信息：
- 花费的时间
- 总运行时间
- 百分比
- 运行计数
如果将鼠标悬停在27行上，它将为您提供以下统计信息：
总结
在进行科学操作时，了解 NumPy 的内部结构至关重要。 效率是关键，因为许多科学计算都是计算和内存密集型的。 因此，如果您的代码编写效率不高，则计算所需的时间将远远超过所需的时间，这将损害您的研发时间。
在本章中，您已经了解了 NumPy 库的一些内部和性能方面，还了解了vprof库，该库可帮助您检查 python 程序的性能。
代码概要分析将极大地帮助您逐行检查程序，并且如您先前所见，查看相同数据的方式也有所不同。 确定了程序中最苛刻的部分后，便可以开始搜索更有效的方法或实现以提高性能并节省更多时间。
在下一章中，我们将概述高性能，低级的数值计算库。 NumPy 可以使用这些实现来获得可观的性能提升。
八、高性能数值计算库概述
在科学计算应用中可以执行许多数值运算，并且未经优化的代码或库实现会导致严重的性能瓶颈。
NumPy 库通过更有效地使用其内存布局来帮助提高 Python 程序的性能。
在实际应用中，最常用的数学分支之一是线性代数。 线性代数用于计算机图形学，密码学，计量经济学，机器学习，深度学习和自然语言处理，仅举几个例子。 具有高效的矩阵和向量运算至关重要。
高性能，低级框架（例如 BLAS，LAPACK 和 ATLAS）是 Netlib 库的一部分，用于密集的线性代数运算；其他框架（例如 Intel MKL）也可以在其中使用您的程序。 这些库在计算中具有很高的性能和准确性。 您可以通过其他高级编程语言（例如 Python 或 C++）调用它们来使用这些库。
当 NumPy 与不同的 BLAS 库链接时，您可以观察到性能差异而无需更改代码，因此了解哪种链接可以更好地提高性能非常重要。
让我们看一下这些低级库。
BLAS 和 LAPACK
BLAS 代表基本线性代数子程序，并且是处理线性代数运算的低级例程的标准。 低级例程包括向量和矩阵加/乘，线性组合等操作。 BLAS 为线性代数运算提供了三种不同的级别：
- BLAS1：标量向量和向量向量运算
- BLAS2：矩阵向量运算
- BLAS3：矩阵矩阵运算
LAPACK 代表线性代数软件包，并包含更高级的操作。 LAPACK 提供了用于矩阵分解（例如 LU，Cholesky 和 QR）以及解决特征值问题的例程。 LAPACK 主要取决于 BLAS 例程。
ATLAS
有许多优化的 BLAS 实现。 ATLAS 代表自动调谐线性代数软件，并且是与平台无关的项目，可以生成优化的 BLAS 实现。
英特尔 MKL
英特尔 MKL 为英特尔处理器优化了 BLAS。 改进了例程和函数，例如 1 级，2 级和 3 级 BLAS，LAPACK 例程，求解器，FFT 函数，其他数学和统计函数。 这些改进的例程和函数得益于共享内存多处理等改进，它们可用于在发行版（如 Anaconda 发行版）中加速科学 python 库（例如 NumPy 和 SciPy）。 如果您查看其发行说明， 可以看到每个发行版都进行了一些重要的改进，例如 LAPACK 函数的性能得到了提高。
OpenBLAS
OpenBLAS 是另一个优化的 BLAS 库，它为不同的配置提供了 BLAS3 级的优化。 作者报告说，与 BLAS 相比，性能增强和改进可与英特尔 MKL 的性能相媲美。
使用 AWS EC2 和底层库配置 NumPy
- 登录到 AWS。 如果您没有帐户，请创建一个：
-
选择 EC2 。
-
单击启动实例：
- 选择
Ubuntu Server 16.04 LTS (HVM), SSD Volume Type - ami-db710fa3：
- 选择
t2.micro实例类型：
- 点击启动：
- 单击启动。
- 选择创建一个新的密钥对：
- 给它命名，然后单击启动实例。 它需要一些时间才能运行：
- 一旦其状态为
running，点击实例 ID ，在这种情况下为i-00ccaeca61a24e042。 然后选择实例并单击Connect：
- 然后它将向您显示以下窗口，其中包含一些有用的信息：
- 打开终端，然后导航到保存所生成密钥的文件夹。 在此示例中，键名称为
aws_oregon。 运行以下命令：
$ chmod 400 aws_oregon.pem
- 然后，在上一个窗口的示例部分中复制该行并运行它：
$ ssh -i "aws_oregon.pem" ubuntu@ec2-34-219-121-1.us-west-2.compute.amazonaws.com
- 在第一个问题的答案中输入
yes将其添加到已知主机中，它将连接到您的实例：
您需要做的第一件事是通过运行以下命令来更新和升级预安装的软件包：
sudo apt-get update
sudo apt-get upgrade
sudo短语为您提供了更新和升级的必要权利，因为软件包的更改可能会对系统产生负面影响，并非所有人都可以对其进行授权。 您可以将apt-get视为 Ubuntu 的软件包管理器。
您可以创建许多虚拟环境，并链接到不同的低级库，但是，每次使用新的低级库配置 NumPy 时，您都将从一个新的预配实例开始。 这将为您提供有关配置过程的想法，以后使您可以轻松地设置虚拟环境。
安装 BLAS 和 LAPACK
为了设置您的开发环境，您需要在运行以下命令后安装所需的软件包，例如编译器，库和其他必要的部分，
$ sudo apt-get update $ sudo apt-get upgrade
对于此配置，很幸运，因为您可以运行以下命令来安装 Python 的 SciPy 软件包，它将安装所有必需的软件包，包括 NumPy，基本线性代数子程序（libblas3）和线性代数软件包（liblapack3）：
$ sudo apt-get install python3-scipy
控制台输出：
- 输入
Y并按Enter继续。 安装完成后，运行以下命令以打开python3解释器：
$ python3
启动 Python 控制台：
- 导入
numpy并使用show_config方法查看 NumPy 的配置：
控制台输出：
- 由于在安装 NumPy 时可以使用 BLAS 和 LAPACK 库，因此它使用它们来构建库。 您可以在
lapack_info和blas_info中看到它们； 否则，您将不会在输出中看到它们，如以下屏幕截图所示：
- 如果您使用的是 macOS，则可以使用 Accelerate/vecLib 框架。 以下命令将输出 macOS 系统的加速器选项：
安装 OpenBLAS
OpenBLAS 的步骤略有不同，如下所示：
- 在先前的配置中运行以下命令：
$ sudo apt-get update $ sudo apt-get upgrade
- 您需要通过运行以下命令来安装
build-essential，其中包括make命令和其他必要的库：
$ sudo apt-get install build-essential libc6 gcc gfortran
- 创建一个名为
openblas_setup.sh的文件，然后粘贴以下内容。 如果您搜索 GitHub，则可以找到不同的设置脚本，并且可以尝试一种满足您需要的脚本：
##!/bin/bash
set -e
pushd /root
git clone https://github.com/xianyi/OpenBLAS.git
pushd /root/OpenBLAS
make clean
make -j4
rm -rf /root/openblas-install
make install PREFIX=/root/openblas-install
popd
ln -sf /root/openblas-install/lib/libopenblas.so /usr/lib/libblas.so
ln -sf /root/openblas-install/lib/libopenblas.so /usr/lib/libblas.so.3
ln -sf /root/openblas-install/lib/libopenblas.so /usr/lib/liblapack.so.3
- 保存此文件并运行以下命令：
$ chmod +777 openblas_setup.sh
$ sudo ./openblas_setup.sh
- 安装完成后，您可以按以下方式安装 numpy 和 scipy：
$ sudo apt-get install python3-pip
$ pip3 install numpy
$ pip3 install scipy
- 现在，您可以像以前一样检查 NumPy 配置：
安装英特尔 MKL
为了针对英特尔 MKL 构建 NumPy 和 SciPy，请按照以下说明进行操作
- 运行以下命令：
$ sudo apt-get update
$ sudo apt-get upgrade
- 您需要安装 Anaconda 发行版，因为 Anaconda 的安装随 Intel MKL 一起提供。 首先使用以下命令下载 Anaconda：
$ wget https://repo.continuum.io/archive/Anaconda3-5.2.0-Linux-x86_64.sh
- 安装完成后，将
cd插入anaconda3/bin并运行python：
$ cd anaconda3/bin
$ ./python
- 您可以像以前一样检查
numpy配置：
安装 ATLAS
为了针对 ATLAS 构建 NumPy，请按照以下说明进行操作
- 运行以下命令：
$ sudo apt-get update
$ sudo apt-get upgrade
- 您需要通过运行以下命令来安装
build-essential，其中包括make命令和其他必要的库：
$ sudo apt-get install build-essential libc6 gcc gfortran
- 然后，您需要安装
atlas：
$ sudo apt-get install libatlas-base-dev
- 您现在可以按以下方式安装
pip和numpy：
$ sudo apt-get install python3-pip
$ pip3 install --no-cache-dir Cython
$ git clone https://github.com/numpy/numpy.git
$ cd numpy
$ cp site.cfg.example site.cfg
$ vi site.cfg
在site.cfg内，您应注释掉地图集行，并将其设置为地图集安装，如下所示：
[atlas]
library_dirs = /usr/local/atlas/lib
include_dirs = /usr/local/atlas/include
然后运行：
$ sudo python3 setup.py install
- 安装完成后，安装
scipy：
$ pip3 install scipy
然后返回到您的主目录，启动python解释器并检查numpy配置，这将为您提供以下输出：
>>> import numpy as np
>>> np.show_config()
atlas_blas_info:
include_dirs = ['/usr/include/atlas']
language = c
library_dirs = ['/usr/lib/atlas-base']
define_macros = [('HAVE_CBLAS', None), ('ATLAS_INFO', '"\\"3.10.2\\""')]
libraries = ['f77blas', 'cblas', 'atlas', 'f77blas', 'cblas']
...
您已经介绍了上述所有低级库的配置。 是时候了解基准测试的计算密集型任务了。
用于基准测试的计算密集型任务
现在，您将能够使用不同的配置（例如是否使用 BLAS/LAPACK，OpenBLAS，ATLAS 和 Intel MKL）对 NumPy 性能进行基准测试。 让我们回顾一下要为基准计算的内容。
矩阵分解
矩阵分解或分解方法涉及计算矩阵的组成部分，以便可以使用它们简化要求更高的矩阵操作。 在实践中，这意味着将您拥有的矩阵分解为多个矩阵，这样，当您计算这些较小矩阵的乘积时，您将获得原始矩阵。 矩阵分解方法的一些示例是奇异值分解（SVD），特征值分解，Cholesky 分解，下上（LU）和 QR 分解。
奇异值分解
SVD 是线性代数中最有用的工具之一。 Beltrami 和 Jordan 发表了有关其使用的几篇论文。 SVD 用于各种应用，例如计算机视觉和信号处理。
如果您具有正方形或矩形矩阵（M），则可以将其分解为矩阵（U），矩阵（V）（计算中使用矩阵转置）和奇异值（d）。
您的最终公式将如下所示：
以下是奇异值分解的说明：
一种简单的数据精简方法是排除该公式中d小到可以忽略不计的部分。
让我们看看如何使用numpy来实现：
import numpy as np
M = np.random.randint(low=0, high=20, size=20).reshape(4,5)
print(M)
## Output
[[18 15 11 13 19]
[ 1 6 8 13 18]
[ 9 7 15 13 10]
[17 15 12 14 12]]
U, d, VT = np.linalg.svd(M)
print("U:n {}n".format(U))
print("d:n {}n".format(d))
print("VT:n {}".format(VT))
U:
[[-0.60773852 -0.22318957 0.5276743 -0.54990921]
[-0.38123886 0.86738201 0.19333528 0.25480749]
[-0.42657252 0.10181457 -0.82343563 -0.36003255]
[-0.55076919 -0.43297652 -0.07832665 0.70925987]]
d:
[56.31276456 13.15721839 8.08763849 2.51997135]
VT:
[[-0.43547429 -0.40223663 -0.40386674 -0.46371223 -0.52002929] [-0.72920427 -0.29835313 0.06197899 0.27638212 0.54682545]
[ 0.11733943 0.26412864 -0.73449806 -0.30022507 0.53557916]
[-0.32795351 0.55511623 -0.3571117 0.56067806 -0.3773643 ]
[-0.39661218 0.60932187 0.40747282 -0.55144258 0.03609177]]
## Setting full_matrices to false gives you reduced form where small values close to zero are excluded
U, d, VT = np.linalg.svd(M, full_matrices=False)
print("U:n {}n".format(U))
print("d:n {}n".format(d))
print("VT:n {}".format(VT))
## Output
U:
[[-0.60773852 -0.22318957 0.5276743 -0.54990921]
[-0.38123886 0.86738201 0.19333528 0.25480749]
[-0.42657252 0.10181457 -0.82343563 -0.36003255]
[-0.55076919 -0.43297652 -0.07832665 0.70925987]]
d:
[56.31276456 13.15721839 8.08763849 2.51997135]
VT:
[[-0.43547429 -0.40223663 -0.40386674 -0.46371223 -0.52002929]
[-0.72920427 -0.29835313 0.06197899 0.27638212 0.54682545]
[ 0.11733943 0.26412864 -0.73449806 -0.30022507 0.53557916]
[-0.32795351 0.55511623 -0.3571117 0.56067806 -0.3773643 ]]
Cholesky 分解
如果您有一个正方形矩阵，也可以应用 Cholesky 分解，将一个矩阵（M）分解为两个三角形矩阵（U和U^T）。 Cholesky 分解可帮助您简化计算复杂性。 可以将其总结为以下公式：
M = U^T U
以下是 Cholesky 分解的说明：
让我们看看如何使用numpy实现它：
from numpy import array
from scipy.linalg import cholesky
M = np.array([[1, 3, 4],
[2, 13, 15],
[5, 31, 33]])
print(M)
## Output
[[ 1 3 4]
[ 2 13 15]
[ 5 31 33]]
L = cholesky(M)
print(L)
## Output
[[1\. 3\. 4\. ]
[0\. 2\. 1.5 ]
[0\. 0\. 3.84057287]]
L.T.dot(L)
## Output
array([[ 1., 3., 4.],
[ 3., 13., 15.],
[ 4., 15., 33.]])
LU 分解
与 Cholesky 分解类似，LU 分解将矩阵（M）分解为下（L）和上（U）三角矩阵。 这也有助于我们简化计算密集型代数。 可以将其总结为以下公式：
M = LU
下面是 LU 分解的说明：
让我们看看如何使用numpy实现它：
from numpy import array
from scipy.linalg import lu
M = np.random.randint(low=0, high=20, size=25).reshape(5,5)
print(M)
## Output
[[18 12 14 15 2]
[ 4 2 12 18 3]
[ 9 19 5 16 8]
[15 19 6 16 11]
[ 1 19 2 18 17]]
P, L, U = lu(M)
print("P:n {}n".format(P))
print("L:n {}n".format(L))
print("U:n {}".format(U))
## Output
P:
[[1\. 0\. 0\. 0\. 0.]
[0\. 0\. 1\. 0\. 0.]
[0\. 0\. 0\. 0\. 1.]
[0\. 0\. 0\. 1\. 0.]
[0\. 1\. 0\. 0\. 0.]]
L:
[[ 1\. 0\. 0\. 0\. 0\. ]
[ 0.05555556 1\. 0\. 0\. 0\. ]
[ 0.22222222 -0.03636364 1\. 0\. 0\. ]
[ 0.83333333 0.49090909 -0.70149254 1\. 0\. ]
[ 0.5 0.70909091 -0.32089552 0.21279832 1\. ]]
U:
[[18\. 12\. 14\. 15\. 2\. ]
[ 0\. 18.33333333 1.22222222 17.16666667 16.88888889]
[ 0\. 0\. 8.93333333 15.29090909 3.16969697]
[ 0\. 0\. 0\. 5.79918589 3.26594301]
[ 0\. 0\. 0\. 0\. -4.65360318]]
P.dot(L).dot(U)
## Output
array([[18., 12., 14., 15., 2.],
[ 4., 2., 12., 18., 3.],
[ 9., 19., 5., 16., 8.],
[15., 19., 6., 16., 11.],
[ 1., 19., 2., 18., 17.]])
特征值分解
特征值分解也是一种适用于平方矩阵的分解技术。 使用特征值分解分解方阵（M）时，将得到三个矩阵。 这些矩阵之一（Q）在列中包含特征向量，另一个矩阵（L）在对角线中包含特征值，最后一个矩阵是特征向量矩阵（Q^(-1)）。
可以将其总结为以下公式：
M = QVQ^(-1)
特征值分解将为您提供矩阵的特征值和特征向量。
下面是特征值分解的说明：
让我们看看如何使用numpy实现它：
from numpy import array
from numpy.linalg import eig
M = np.random.randint(low=0, high=20, size=25).reshape(5,5)
print(M)
## Output
[[13 9 5 0 12]
[13 6 11 8 15]
[16 17 15 12 1]
[17 8 5 7 5]
[10 6 18 5 19]]
V, Q = eig(M)
print("Eigenvalues:n {}n".format(V))
print("Eigenvectors:n {}".format(Q))
## Output
Eigenvalues:
[50.79415691 +0.j 5.76076687+11.52079216j
5.76076687-11.52079216j -1.15784533 +3.28961651j
-1.15784533 -3.28961651j]
Eigenvectors:
[[ 0.34875973+0.j -0.36831427+0.21725348j -0.36831427-0.21725348j
-0.40737336-0.19752276j -0.40737336+0.19752276j]
[ 0.46629571+0.j -0.08027011-0.03330739j -0.08027011+0.03330739j
0.58904402+0.j 0.58904402-0.j ]
[ 0.50628483+0.j 0.62334823+0.j 0.62334823-0.j
-0.27738359-0.22063552j -0.27738359+0.22063552j]
[ 0.33975886+0.j 0.14035596+0.39427693j 0.14035596-0.39427693j
0.125282 +0.46663129j 0.125282 -0.46663129j]
[ 0.53774952+0.j -0.18591079-0.45968785j -0.18591079+0.45968785j
0.20856874+0.21329768j 0.20856874-0.21329768j]]
from numpy import diag
from numpy import dot
from numpy.linalg import inv
Q.dot(diag(V)).dot(inv(Q))
## Output
array([[1.30000000e+01-2.88657986e-15j, 9.00000000e+00-2.33146835e-15j,
5.00000000e+00+2.38697950e-15j, 1.17683641e-14+1.77635684e-15j,
1.20000000e+01-4.99600361e-16j],
[1.30000000e+01-4.32986980e-15j, 6.00000000e+00-3.99680289e-15j,
1.10000000e+01+3.38618023e-15j, 8.00000000e+00+1.72084569e-15j,
1.50000000e+01-2.77555756e-16j],
[1.60000000e+01-7.21644966e-15j, 1.70000000e+01-6.66133815e-15j,
1.50000000e+01+5.71764858e-15j, 1.20000000e+01+2.99760217e-15j,
1.00000000e+00-6.66133815e-16j],
[1.70000000e+01-5.27355937e-15j, 8.00000000e+00-3.10862447e-15j,
5.00000000e+00+4.27435864e-15j, 7.00000000e+00+2.22044605e-15j,
5.00000000e+00-1.22124533e-15j],
[1.00000000e+01-3.60822483e-15j, 6.00000000e+00-4.21884749e-15j,
1.80000000e+01+2.27595720e-15j, 5.00000000e+00+1.55431223e-15j,
1.90000000e+01+3.88578059e-16j]])
QR 分解
您可以通过应用 QR 分解将正方形或矩形矩阵（M）分解为正交矩阵（Q）和上三角矩阵（R）。 可以用以下公式表示：
M = QR
以下是 QR 分解的说明：
让我们看看如何使用numpy实现它：
from numpy import array
from numpy.linalg import qr
M = np.random.randint(low=0, high=20, size=20).reshape(4,5)
print(M)
## Output
[[14 6 0 19 3]
[ 9 6 17 8 8]
[ 4 13 17 4 4]
[ 0 0 2 7 11]]
Q, R = qr(M, 'complete')
print("Q:n {}n".format(Q))
print("R:n {}".format(R))
## Output
Q:
[[-0.81788873 0.28364908 -0.49345895 0.08425845]
[-0.52578561 -0.01509441 0.83834961 -0.14314877]
[-0.2336825 -0.95880935 -0.15918031 0.02718015]
[-0\. -0\. 0.16831464 0.98573332]]
R:
[[-17.11724277 -11.09991852 -12.91095786 -20.68090082 -7.59468109]
[ 0\. -10.85319349 -16.5563638 1.43333978 -3.10504542]
[ 0\. 0\. 11.88250752 -2.12744187 6.4411599 ]
[ 0\. 0\. 0\. 7.4645743 10.05937231]]
array([[1.40000000e+01, 6.00000000e+00, 1.77635684e-15, 1.90000000e+01,
3.00000000e+00],
[9.00000000e+00, 6.00000000e+00, 1.70000000e+01, 8.00000000e+00,
8.00000000e+00],
[4.00000000e+00, 1.30000000e+01, 1.70000000e+01, 4.00000000e+00,
4.00000000e+00],
[0.00000000e+00, 0.00000000e+00, 2.00000000e+00, 7.00000000e+00,
1.10000000e+01]])
处理稀疏线性系统
您将不会总是使用密集矩阵，并且当您需要使用稀疏矩阵时，有些库将帮助您优化稀疏矩阵运算。 即使这些可能没有 Python API，您仍可能需要通过使用其他编程语言来使用它们，例如 C 和 C++：
- Hypre：包含预处理器和求解器，以利用并行实现来处理稀疏线性方程组。
- SuperLU：处理大型，稀疏，不对称的线性方程组。
- UMFPACK：解决稀疏线性方程组。
- CUSP：带有并行实现的稀疏线性代数和图形计算的开源库。 通过使用 CUSP，您可以访问 NVIDIA GPU 提供的计算资源。
- cuSPARSE：包含用于处理稀疏矩阵的线性代数子例程。 与 CUSP 一样，您可以访问 Nvidia GPU 提供的计算资源。
总结
在本章中，您探索了可以与 NumPy 配对的各种低级库及其配置。 我们特意运行了 EC2 条款，以便您熟悉基本的 Linux 命令行操作。 您还研究了各种计算密集型，数值，线性代数运算，这些运算将在下一章中用作基准测试不同的配置。
在下一章中，我们将创建一个基准 python 脚本，以在每种配置上运行。 您将能够查看不同线性代数运算和不同矩阵大小的性能指标
