原文：Mastering Numerical Computing With NumPy

协议：CC BY-NC-SA 4.0

译者：飞龙

六、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) 

您可以从前面的代码块中获得以下图表：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//545138a3-1597-4e73-ac37-c900364a3514.png

使用以下代码将图像数组转换为二维数据集：

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)) 

您从前面的代码中获得以下输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//2b717041-1181-43d5-9a52-88acdfba12d0.png

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() 

我们从前面代码的输出中获得以下图表，如下图所示：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//689cca38-243b-492a-886a-0f3b418eeed8.png

您还可以查看平均房间数与房价之间的关系。 以下代码块打印出性能指标：

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() 

我们从前面的代码中获得以下输出，如下图所示：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//28064f39-39a5-4056-be40-7912e65b6dcf.png

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) 

我们从前面的代码中获得以下输出，如下表所示：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//bde1f31c-2242-48f3-9595-9b5573236667.png

这段代码向您展示了如何在数据帧中创建目标列：

df['Target'] = diabetes.target 
df.head(10) 

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//05017bb6-90a0-4304-99bc-b92f2a81d26c.png

Pandas 帮助我们轻松地处理表格数据，并通过各种辅助方法和可视化支持我们的分析。 看一下代码：

## Descriptive statistics 
df.describe() 

我们从前面的代码中获得以下输出，如下表所示：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//c4ae24f8-f89e-42c4-b444-731c70aaaa26.png

通过使用以下代码行，看看目标的分布方式：

plt.hist(df['Target']) 

下图显示了上一行的输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//cc3d1506-c053-4b96-9664-5f6996dcc72b.png

您可以看到目标变量向右倾斜。 看一下这段代码：

## 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 

前一个代码块的输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//8c302060-2632-4875-ba97-297f22708c38.png

Feature distributions

您可以检查其中一些特征的分布，并确定其中哪些看起来相似。 对于此示例，特征s1，s2和s6似乎具有相似的分布，如从此代码中可以看到的：

## corr method will give you the correlation between features 
df[numeric_cols].corr() 

下图显示了上一行的输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//432ec7dc-28c7-42b2-aec2-4e1cb680b93c.png

您可以使用heatmap更好地表示这种关系，如下所示：

plt.figure(figsize=(15, 15)) 
sns.heatmap(df[numeric_cols].corr(), annot=True) 

下图是由前面的代码块生成的热图：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//3052dc7f-46c9-40d8-b086-1d681af5fbb6.png

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) 

此图显示了过滤后的相关性：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//9a1ecb64-ffa6-419d-bbde-3dc360fddc4a.png

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) 

该图显示了前面代码的以下输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//0667393f-c07d-4474-ab67-85ba7c7b7aa1.png

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()的输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//886d9310-a3b9-4f9d-b4d8-3ee461df7f72.png

让我们通过导入以下设置来自定义绘图：

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()

该图显示了来自先前设置的图：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//329ac184-fb7f-46c2-aa1f-9c910e9b692d.png

您可以使用以下代码来计算每日变化：

msft['Daily Pct. Change'] = (msft['Adj. Close'] - msft['Adj. Open']) / msft['Adj. Open'] 

msft.tail(10)

下图显示了msft.tail(10)的输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//1f3615cd-31a2-4102-bb3d-ede8046d35cd.png

尝试使用此每日收益的直方图：

plt.figure(figsize=(22, 8)) 
plt.hist(msft['Daily Pct. Change'], bins=100)

如下图所示，它将为您提供以下图表：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//4a4095c3-b531-4704-b8ed-2bfc803dcee4.png

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)

我们从前面的代码中获得以下图表，如下图所示：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//d68b592e-855e-48d0-846c-43277a2935c7.png

您可以创建开盘-最高-最低-收盘价（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)

该图显示了先前代码的图：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//a5cf9226-e77d-4bb5-908b-7c11c7241c00.png

您可以通过在图表上选择自定义范围来检查特定区域。 看一下这个图：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//59f2162c-bdf9-426b-a3d2-216a3bc94d4b.png

同样，您可以使用以下代码创建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)

下图显示了此代码的输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//9d1f6097-f58a-4ff8-8e1f-db34ad7f3ae4.png

您也可以在Candlestick图表中选择特定范围。 看一下这个图：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//39a01243-e40b-4aff-a059-db2b42e33b07.png

分布图如下：

import plotly.figure_factory as ff 

fig = ff.create_distplot([msft['Daily Pct. Change'].values], ['MSFT Daily Returns'], show_hist=False) 

plot(fig)

下图显示了前面代码的输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//d49a8f10-1559-47d2-860b-83faf46b9716.png

我们可以创建三个移动平均线，如下所示：

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)的输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//f3050952-395c-49ad-8076-64fa13ba7813.png

根据切片的数据，将包括最近 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)

下图显示了前面代码块中的图：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//242fa07e-9351-4d35-8078-8d80075930ca.png

移动平均线用于监视金融市场的趋势。 在此示例中，有三个移动平均线，每个移动线均具有不同的周期。 您可以设置分析天数，以进行短期，中期和长期趋势监视。

当您开始使用财务时间序列时，您会很快意识到您需要基于不同期间的汇总，并且在pandas中创建这些汇总非常容易。 以下代码段将通过计算平均值每月汇总记录：

msft_monthly = msft.resample('M').mean() 

msft_monthly.tail(10)

下图显示了msft_monthly.tail(10)的输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//e54a85cb-898d-4b7a-a9eb-996e8360618d.png

这是一个简单的时间序列图：

data = [go.Scatter(x=msft_monthly[-24:].index, y = msft_monthly[-24:]['Adj. Close'])] 

plot(data)

如下图所示，这将为您提供以下图表：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//8d35cdbf-b546-43ce-bfa4-5b2d2cc1f17c.png

在检查要素之间的关系时，可以使用我们在前面的示例中已经看到的相关矩阵。 在时间序列中，从业者对自相关感兴趣，自相关显示了时间序列与其自身的相关性。 例如，理想情况下，您希望时间序列中出现周期性的峰值，以显示您的季节性。 通过使用以下代码，让我们看看每日百分比变化是否有任何明显的峰值：

plt.figure(figsize=(22, 14)) pd.plotting.autocorrelation_plot(msft_monthly['Daily Pct. Change'])

我们从前面的代码中获得以下图表，如下图所示：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//4f88b706-cb1a-465a-95b5-cfe7c8802673.png

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]}

下图显示了数据类型树：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//490b0e1b-692d-4f2a-855d-f3482dd0da00.png

您可以通过调用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]

将会选择：

1. 首先，第一维的所有项目
2. 然后，第二维的每两个项目
3. 然后，第三维的每两个项目

它将具有以下数组：

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...

一个新的选项卡将在您的浏览器中打开，并显示以下输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//30c1b63e-e750-467f-89f8-690435d78992.png

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...

新标签页将在您的浏览器中打开，并显示以下输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//c47f02ee-21bf-484c-abd0-727be115ea3f.png

内存使用情况

在左侧，您可以看到内存中的对象，并且图表显示了随着执行的行数增加，内存使用量（以兆字节为单位）。 如果将鼠标悬停在图表上，则每个点将具有以下信息：

• 执行的行
• 行号
• 函数名称
• 文件名称
• 内存使用情况

例如，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...

“新”标签将在您的浏览器中打开，并显示以下输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//0f7b1b59-ce7f-4ec5-930d-af4f4923ff5f.png

Heatmap for the lines executed

在左侧，有一个已检查模块的列表，在这种情况下，只有一个文件要检查。

右边是程序中每行的热图。 如果将鼠标悬停在行上，它将为您提供以下信息：

• 花费的时间
• 总运行时间
• 百分比
• 运行计数

如果将鼠标悬停在27行上，它将为您提供以下统计信息：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//2a01ed56-5ab9-429c-923f-b9880202e834.png

总结

在进行科学操作时，了解 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

1. 登录到 AWS。 如果您没有帐户，请创建一个：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//fa6d4e79-b597-418d-8ecc-c8c3e8dbec9d.png

2. 选择 EC2 。

3. 单击启动实例：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//dfeff87d-2771-46bd-bd75-c53c73d0ce97.png

4. 选择Ubuntu Server 16.04 LTS (HVM), SSD Volume Type - ami-db710fa3：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//b4cda929-039a-4c1c-973e-666143148257.png

5. 选择t2.micro实例类型：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//33e941fa-99d4-43e7-b93b-373f7cb00455.png

6. 点击启动：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//151e7c10-aef3-40da-923a-b7f7171dbbf5.png

7. 单击启动。
8. 选择创建一个新的密钥对：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//dc5df025-8ebf-4109-a0d9-5fc104043bc6.png

9. 给它命名，然后单击启动实例。 它需要一些时间才能运行：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//41e72212-d346-45d0-bf0d-6d42645393da.png

10. 一旦其状态为running，点击实例 ID ，在这种情况下为i-00ccaeca61a24e042。 然后选择实例并单击Connect：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//d207b825-c27d-4e72-817c-e7256a846d7e.png

11. 然后它将向您显示以下窗口，其中包含一些有用的信息：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//70613b58-c8f1-4488-a96e-a852ccf32900.png

12. 打开终端，然后导航到保存所生成密钥的文件夹。 在此示例中，键名称为aws_oregon。 运行以下命令：

$ chmod 400 aws_oregon.pem

13. 然后，在上一个窗口的示例部分中复制该行并运行它：

$ ssh -i "aws_oregon.pem" ubuntu@ec2-34-219-121-1.us-west-2.compute.amazonaws.com

14. 在第一个问题的答案中输入yes将其添加到已知主机中，它将连接到您的实例：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//2f4724ca-92f2-4ea2-badd-edbd3a32513a.png

您需要做的第一件事是通过运行以下命令来更新和升级预安装的软件包：

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

控制台输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//0105fda2-1f28-41df-9f3b-3ae8857ce411.png

1. 输入Y并按Enter继续。 安装完成后，运行以下命令以打开python3解释器：

$ python3

启动 Python 控制台：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//782fd1c3-0c4d-40da-9c2f-456bf1011bb2.png

2. 导入numpy并使用show_config方法查看 NumPy 的配置：

控制台输出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//961410ae-8b32-4f4c-894c-736a7e2a020b.png

3. 由于在安装 NumPy 时可以使用 BLAS 和 LAPACK 库，因此它使用它们来构建库。 您可以在lapack_info和blas_info中看到它们； 否则，您将不会在输出中看到它们，如以下屏幕截图所示：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//f5e7565b-79d8-4432-b04a-7d19c130a123.png

4. 如果您使用的是 macOS，则可以使用 Accelerate/vecLib 框架。 以下命令将输出 macOS 系统的加速器选项：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//f5f87ec3-4393-48fb-8491-3953cbb2a9ef.png

安装 OpenBLAS

OpenBLAS 的步骤略有不同，如下所示：

1. 在先前的配置中运行以下命令：

$ sudo apt-get update $ sudo apt-get upgrade

2. 您需要通过运行以下命令来安装build-essential，其中包括make命令和其他必要的库：

$ sudo apt-get install build-essential libc6 gcc gfortran

3. 创建一个名为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

4. 保存此文件并运行以下命令：

$ chmod +777 openblas_setup.sh
$ sudo ./openblas_setup.sh

5. 安装完成后，您可以按以下方式安装 numpy 和 scipy：

$ sudo apt-get install python3-pip
$ pip3 install numpy
$ pip3 install scipy

6. 现在，您可以像以前一样检查 NumPy 配置：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//efaae3d9-8ade-4949-b31b-895030db9223.png

安装英特尔 MKL

为了针对英特尔 MKL 构建 NumPy 和 SciPy，请按照以下说明进行操作

1. 运行以下命令：

$ sudo apt-get update
$ sudo apt-get upgrade

2. 您需要安装 Anaconda 发行版，因为 Anaconda 的安装随 Intel MKL 一起提供。 首先使用以下命令下载 Anaconda：

$ wget https://repo.continuum.io/archive/Anaconda3-5.2.0-Linux-x86_64.sh

3. 安装完成后，将cd插入anaconda3/bin并运行python：

$ cd anaconda3/bin
$ ./python

4. 您可以像以前一样检查numpy配置：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//c8d03c39-881f-422a-b592-34a20dd532ea.png

安装 ATLAS

为了针对 ATLAS 构建 NumPy，请按照以下说明进行操作

1. 运行以下命令：

$ sudo apt-get update
$ sudo apt-get upgrade

2. 您需要通过运行以下命令来安装build-essential，其中包括make命令和其他必要的库：

$ sudo apt-get install build-essential libc6 gcc gfortran

3. 然后，您需要安装atlas：

$ sudo apt-get install libatlas-base-dev

4. 您现在可以按以下方式安装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

5. 安装完成后，安装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）。

您的最终公式将如下所示：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//5d580481-4701-4532-a958-660363f75c69.png

以下是奇异值分解的说明：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//7afa6d97-21e5-40bf-8138-43a5c61c5f69.png

一种简单的数据精简方法是排除该公式中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 分解的说明：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//6344c2d0-3b1f-442b-bdd4-302236b86deb.png

让我们看看如何使用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 分解的说明：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//9a2c2918-c0a7-4f19-abce-bbe8a14b3fa1.png

让我们看看如何使用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)

特征值分解将为您提供矩阵的特征值和特征向量。

下面是特征值分解的说明：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//cb9a2709-fbfc-4fc0-be93-10fcbf54b757.png

让我们看看如何使用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 分解的说明：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//ac3190aa-d5a4-449c-8b17-43ab32e841da.png

让我们看看如何使用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 脚本，以在每种配置上运行。 您将能够查看不同线性代数运算和不同矩阵大小的性能指标

九、性能基准

在本章中，您将研究上一章介绍的不同配置的性能统计信息。 当然，当前设置无法为您提供最准确的环境，因为您无法控制 EC2 实例，但是它将使您了解自己环境中所需的设置。

我们将涵盖以下主题：

• 基准的重要性
• BLAS，LAPACK，OpenBLAS，ATLAS 和 Intel MKL 的性能
• 最终结果

为什么我们需要基准？

随着编程技巧的提高，您将开始实施更高效的程序。 您将搜索数十个代码存储库，以了解其他人如何解决类似的问题，并且您会发现那些令您赞叹不已的稀有宝石。

在编写更好的软件和实施系统的整个过程中，您将需要测量和跟踪改进速度的方法。 通常，您将以起点为基准，并查看所做的改进将如何累加性能指标。

设置基准后，您将对几种不同的实现进行基准测试，并将有机会根据您选择的性能指标进行比较。 您可以选择各种指标，并且需要事先决定。

这些基准的性能指标将保持相当简单，并且仅使用所花费的时间指标。 您将使用不同的配置多次执行相同的操作，并首先计算平均花费的时间。 计算平均值的公式为：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//c928d7f5-e3dd-4030-85f1-289acd24f6d5.png

这是一个计算均值的好公式； 在我们的示例中，公式解释如下：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//123d5f4e-e787-49e3-844f-23b6dbe90944.png

基线将基于此公式创建。 第一组计算如下：

的加法和乘法：

• 向量向量
• 向量矩阵
• 矩阵矩阵

通常，您将运行这些计算给定次数并计算平均值。

以下代码段向您展示了一个自定义函数，而不是 Python 中可用的通用计时器。 使用自定义函数的原因是，您以后可以将其与其他统计函数一起扩展，并通过适当的日志记录更好地查看详细信息。 函数将在计算开始之前输出有用的信息，并在迭代完成之后输出结果。

import inspect
import time
from datetime import datetime

def timer(*args, operation, n):
 """
 Returns average time spent 
 for given operation and arguments.

 Parameters
 ----------
 *args: list (of numpy.ndarray, numpy.matrixlib.defmatrix.matrix or both)
 one or more numpy vectors or matrices
 operation: function
 numpy or scipy operation to be applied to given arguments
 n: int 
 number of iterations to apply given operation
 Returns
 -------
 avg_time_spent: double
 Average time spent to apply given operation
 std_time_spent: double
 Standard deviation of time spent to apply given operation

 Examples
 --------

 >>> import numpy as np
 >>> vec1 = np.array(np.random.rand(1000))
 >>> vec2 = np.array(np.random.rand(1000))
 >>> args = (vec1, vec2)

 >>> timer(*args, operation=np.dot, n=1000000)
 8.942582607269287e-07
 """

 # Following list will hold the
 # time spent value for each iteration
 time_spent = []

 # Configuration info
 print("""
 -------------------------------------------

 ### {} Operation ###

 Arguments Info
 --------------
 args[0] Dimension: {},
 args[0] Shape: {},
 args[0] Length: {}
 """.format(operation.__name__,
 args[0].ndim,
 args[0].shape,
 len(args[0])))

 # If *args length is greater than 1, 
 # print out the info for second argument
 args_len = 0
 for i, arg in enumerate(args):
     args_len += 1

 if args_len > 1:
     print("""
     args[1] Dimension: {},
     args[1] Shape: {},
     args[1] Length: {}
     """.format(args[1].ndim,
         args[1].shape,
         len(args[1])))

 print("""
 Operation Info
 --------------
 Name: {},
 Docstring: {}

 Iterations Info
 ---------------
 # of iterations: {}""".format(
 operation.__name__,
 operation.__doc__[:100] + 
 "... For more info type 'operation?'",
 n))

 print("""
 -> Starting {} of iterations at: {}""".format(n, datetime.now()))

 if args_len > 1:
     for i in range(n):
         start = time.time()
         operation(args[0], args[1])
         time_spent.append(time.time()-start)
 else:
     for i in range(n):
         start = time.time()
         operation(args[0])
         time_spent.append(time.time()-start)

 avg_time_spent = np.sum(time_spent) / n
 print("""
 -> Average time spent: {} seconds,
 -------------------------------------------
 """.format(avg_time_spent))

 return avg_time_spent

当此函数中包含Docstring时，可以显示它以查看函数参数，返回的内容以及用法示例：

print(timer.__doc__)

这将生成以下输出：

    Returns average time spent 
    for given operation and arguments.

    Parameters
    ----------
        *args: list (of numpy.ndarray, numpy.matrixlib.defmatrix.matrix or both)
            one or more numpy vectors or matrices
        operation: function
            numpy or scipy operation to be applied to given arguments
        n: int 
            number of iterations to apply given operation

    Returns
    -------
        avg_time_spent: double
            Average time spent to apply given operation

    Examples
    --------
    >>> import numpy as np

    >>> vec1 = np.array(np.random.rand(1000))
    >>> vec2 = np.array(np.random.rand(1000))

    >>> args = [vec1, vec2]

    >>> timer(*args, operation=np.dot, n=1000000)
    8.942582607269287e-07

让我们开始测量两个向量的点积所花费的平均时间。 以下代码块定义向量，并创建要输入到计时器函数中的参数：

import numpy as np
vec1 = np.array(np.random.rand(1000))
vec2 = np.array(np.random.rand(1000))
args = [vec1, vec2]

您现在可以按以下方式调用计时器函数：

timer(*args, operation=np.dot, n=1000000)
    -------------------------------------------
    ### dot Operation ###
    Arguments Info
    --------------
    args[0] Dimension: 1,
    args[0] Shape: (1000,),
    args[0] Length: 1000
    args[1] Dimension: 1,
    args[1] Shape: (1000,),
    args[1] Length: 1000
    Operation Info
    --------------
    Name: dot,
    Docstring: dot(a, b, out=None)
    Dot product of two arrays. Specifically,
    - If both `a` and `b` are 1-D... For more info type 'operation?'
    Iterations Info
    ---------------
    # of iterations: 1000000
    -> Starting 1000000 of iterations at: 2018-06-09 21:02:51.711211 
    -> Average time spent: 1.0054986476898194e-06 seconds,
    -------------------------------------------
1.0054986476898194e-06

我们的向量乘积平均需要 1 微秒。 让我们看看如何通过添加其他指标来改进此计算。 您可以轻松添加的另一个指标是标准差，如下公式所示：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//ac8c49f8-3be4-4b3b-82bc-5bc2a6eacd0a.png

您熟悉上图中的公式术语。 标准差只是告诉您所报告指标的可变性，这是在我们的示例中花费的平均时间。

通过计算std_time_spent，打印其值并返回以下内容来扩展计时器函数：

avg_time_spent = np.sum(time_spent) / n
std_time_spent = np.std(time_spent)
print("""
-> Average time spent: {} seconds,
-> Std. deviation time spent: {} seconds
""".format(avg_time_spent, std_time_spent))
return avg_time_spent, std_time_spent

您还可以如下更新Docstring：

Returns
-------
avg_time_spent: double
Average time spent to apply given operation
std_time_spent: double
Standard deviation of time spent to apply given operation.

您可以重新定义时间函数，然后再次运行先前的计算，如下所示：

timer(*args, operation=np.dot, n=1000000)

您将获得以下输出（为简洁起见，仅显示最后一部分）以及其他信息：

-> Starting {} of iterations at: {}".format(n, datetime.now())
-> Average time spent: 1.0006928443908692e-06 seconds,
-> Std. deviation time spent: 1.2182541822530471e-06 seconds
(1.0006928443908692e-06, 1.2182541822530471e-06)

大！ 您还添加了哪些其他指标？ 如何添加置信区间？ 该部分将留给您锻炼，但是对您来说应该很容易！

让我们继续向量矩阵乘积：

mat1 = np.random.rand(1000,1000)
args = [vec1, mat1]
timer(*args, operation=np.dot, n=1000000)

这将为您提供以下输出：

    Arguments Info
    --------------
    args[0] Dimension: 1,
    args[0] Shape: (1000,),
    args[0] Length: 1000
    args[1] Dimension: 2,
    args[1] Shape: (1000, 1000),
    args[1] Length: 1000
    Operation Info
    --------------
    Name: dot,
    Docstring: dot(a, b, out=None)
    Dot product of two arrays. Specifically,
    - If both `a` and `b` are 1-D... For more info type 'operation?'
    Iterations Info
    ---------------
    # of iterations: 1000000
    -> Starting 1000000 of iterations at: 2018-06-09 19:13:07.013949
    -> Average time spent: 0.00020063393139839174 seconds,
    -> Std. deviation time spent: 9.579314466482879e-05 seconds
 (0.00020063393139839174, 9.579314466482879e-05)

最后，矩阵矩阵乘法如下：

mat1 = np.random.rand(100,100)
mat2 = np.random.rand(100,100)
args = [mat1, mat2]
timer(*args, operation=np.dot, n=1000000)

这将为您提供类似于先前输出的输出。

现在，我们或多或少有了一个想法，即如何挑战如何在计算机上执行这些任务。 基准函数列表已完成，在上一章中您看到了将点积添加到矩阵分解中的信息。

您将要做的是创建一个包含这些计算和统计信息的 Python 脚本文件。 然后，您将使用在 AWS 上设置的不同配置运行此文件。

让我们看一下linalg_benchmark.py，您可以在这个页面 中找到它。

以下代码块向您展示了linalg_benchmark.py脚本的重要部分，该脚本将用于测试您先前在 AWS 上设置的不同配置：

## Seed for reproducibility
np.random.seed(8053)
dim = 100
n = 10000
v1, v2 = np.array(rand(dim)), np.array(rand(dim))
m1, m2 = rand(dim, dim), rand(dim, dim)
## Vector - Vector Product
args = [v1, v2]
timer(*args, operation=np.dot, n=n)
## Vector - Matrix Product
args = [v1, m1]
timer(*args, operation=np.dot, n=n)
## Matrix - Matrix Product
args = [m1, m2]
timer(*args, operation=np.dot, n=n)
## Singular-value Decomposition
args = [m1]
timer(*args, operation=np.linalg.svd, n=n)
## LU Decomposition
args = [m1]
timer(*args, operation=lu, n=n)
## QR Decomposition
args = [m1]
timer(*args, operation=qr, n=n)
## Cholesky Decomposition
M = np.array([[1, 3, 4],
[2, 13, 15],
[5, 31, 33]])
args = [M]
timer(*args, operation=cholesky, n=n)
## Eigenvalue Decomposition
args = [m1]
timer(*args, operation=eig, n=n)
print("""
NumPy Configuration:
--------------------
""")
np.__config__.show()

将有两个单独的运行：

• 1 ^st 与dim = 100一起运行
• 2 ^和与dim = 500一起运行

让我们看一下结果。

准备性能基准

对于每个实例和配置，导航到您的Home目录并创建一个名为py_scripts的文件夹：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//71b5bdc9-fabb-4070-9985-89ceab33443f.png

使用以下命令创建名为linalg_benchmark.py的文件并粘贴内容：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//9aa9878f-ace1-4156-9c40-fb98f19ddb67.png

粘贴内容后，键入:，然后键入wq!和Enter保存并退出：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//fef3b1d6-8f81-4b71-8f58-b915c9ffb4f7.png

现在，您可以使用以下命令运行该文件：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//aec20b2e-8717-429d-928b-2d3132bc3722.png

对于 Anaconda 分发，您将使用以下命令运行脚本：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//7ded6a60-0e2b-4702-9223-56f6258b6d60.png

BLAS 和 LAPACK 的性能

在这里，您将使用 BLAS 和 LAPACK 运行linalg_benchmark.py脚本。 连接到具有此配置的t2.micro实例，然后如上一节中所示运行脚本。

以下是dim = 100的运行结果：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//da9e67f5-f8b7-4d54-bc7e-c7c79c8cdee4.png

以下是dim = 500的运行结果：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//0c368f9a-ec08-4b18-b19b-c624c86bcc8c.png

OpenBLAS 的性能

在这里，您将使用 OpenBLAS 运行linalg_benchmark.py脚本。 连接到具有此配置的t2.micro实例，然后运行上一节中显示的脚本。

以下是dim = 100的运行结果：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//90294475-3865-488a-ad81-774c87fda1ee.png

以下是dim = 500的运行结果：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//c09d200e-7fec-48f8-a80f-29c55be47084.png

ATLAS 的性能

在这里，您将使用 ATLAS 运行linalg_benchmark.py脚本。 连接到具有此配置的t2.micro实例，如上一节中所示运行脚本。

以下是dim = 100的运行结果：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//8e4ceff4-cabb-4c65-b790-c9b367d4d195.png

以下是dim = 500的运行结果：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//9097e7ac-a57d-4a27-b5e6-4650b9e84ad6.png

英特尔 MKL 的性能

在这里，您将使用英特尔 MKL 运行linalg_benchmark.py脚本。 连接到具有此配置的t2.micro实例，然后运行上一节中显示的脚本。

以下是dim = 100的运行结果：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//a372c4c3-7c0b-4cd8-98b5-a4891534fec2.png

以下是dim = 500的运行结果：

https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-num-comp-numpy/img//a922450d-3125-4e60-8c76-fa64a1ac4986.png

结果

当然，t2.micro实例相当薄弱，您应该更多地了解 Amazon 如何为 EC2 实例提供这种计算能力。 您可以在这个页面 上阅读有关它们的更多信息。

如果您使用功能更强大的计算机并具有更多的内核，则不同配置之间的性能差异将更加明显。

说到结果，毫不奇怪，默认安装的 BLAS 和 LAPACK 为我们提供了基准性能，而经过优化的版本（如 OpenBLAS，ATLAS 和 Intel MKL）提供了更好的性能。

正如您已经指出的那样，您没有在 Python 脚本中更改任何代码行，而仅通过将 NumPy 库与不同的加速器链接起来，便获得了巨大的性能提升。

如果您将更深入地研究这些低级库以了解提供了哪些特定的例程和函数，则将更好地了解程序的哪些部分将从这些实现中受益。

当然，起初您可能还不了解许多其他细节。 可能是您使用的函数未使用低级库或未并行化操作的情况。 在某些情况下，多线程会或不会有所帮助。 知识和经验最终取决于您的实验，并且您将从自己的经验中学到东西，因此您将更加精通各种应用。

许多研究人员发表了实验的设计和结果。 Google 的快速搜索将为您提供大量资源，以阅读和了解这些库在不同硬件和软件配置下的性能。

总结

在本章中，您探讨了执行计算密集型线性代数运算时不同配置的性能。

基准测试是一项严肃的工作，您至少现在已经具备运行基准测试的基本技能。 您在本章中学习的材料远未完成，但是它为您提供了从哪里开始的想法，并且您肯定可以在许多方面进行改进。

您可以看到的一件事是，逐渐增加向量和矩阵的大小时性能指标的行为。 理想情况下，您将需要功能更强大的硬件，但是t2.micro实例在大多数情况下是免费的，或者提供的价格非常便宜。

由于您将需要处理更多的计算密集型工作负载，因此重要的是要了解您的选择以及哪种选择将为您带来最佳性能。 您可以运行这些简单的实验，至少对性能有所了解，这将为您带来很多帮助，并节省时间和金钱。

如果您走了这么远，恭喜！ 我们认为，遍历所有章节并学习相关材料可以提高您在 Python 科学堆栈方面的技能。

我们希望您喜欢阅读这本书，并感谢您的宝贵时间。