Rust机器学习之ndarray

众所周知，Python之所以能成为机器学习的首选语言，与其丰富易用的库有很大关系。某种程度上可以说是诸如numpy、pandas、scikit-learn、matplotlib、pytorch、networks...等一系列科学计算和机器学习库成就了Python今天编程语言霸主的地位。基本上今天的机器学习任务主要就建立在上面列举的这6个库上。这6个库在Rust上都有对应的替代方案。从今天开始，我将撰写系列教程，带大家一起学习如何使用Rust及其库替代Python来更好得完成机器学习任务。

本文是“Rust替代Python进行机器学习”系列文章的第一篇，其他教程请参考下面表格目录：

数据和算法工程师偏爱Jupyter，为了跟Python保持一致的工作环境，文章中的示例都运行在Jupyter上。因此需要各位搭建Rust交互式编程环境（让Rust作为Jupyter的内核运行在Jupyter上），相关教程请参考 Rust交互式编程环境搭建

什么是ndarray？

ndarray是Rust生态中用于处理数组的库。它包含了所有常用的数组操作。简单地说ndarray相当于Rust的numpy。

除了数组操作，ndarray还通过派生包提供其他丰富功能，比如

• ndarray-linalg 用于线性代数运算；
• ndarray-rand 用于产生随机数；
• ndarray-stats 用于统计计算；

可以说ndaary不但包含了numpy的功能，还包含了部分scipy的功能。

值得一提的是，ndarray还很好地支持很多外部特性。比如可以支持rayon做并行计算，支持BLAS进行底层运算加速。因此ndarray的性能非常彪悍。

ndarray 的 BLAS是通过blas-src实现的，blas-src为ndarray提供可选的BLAS源，目前支持

• accelerate - 苹果的Accelerate 框架（仅支持Mac系统）
• blis - BLIS是一个类BLAS高性能线性代数库
• intel-mkl - Intel-MKL是⼀套经过高度优化和广泛线程化的数学库
• netlib - Netlib是一系列数学工具的集合，其中包含BLAS和LAPACK实现
• openblas - OpenBLAS是一套高度优化的开源BLAS库

为什么需要ndarray？

Rust本身已经自带了数组（或列表）类型，也有vector容器。得益于Rust自身的强大，这些自带的数据类型或容器比很多其他编程语言都要快，甚至快一个数量级。那么我们为什么还需要ndarray呢？

首先，ndarray是专门为处理n维数组（矩阵）而设计的，里面包含了很多数学运算，比如矩阵相乘、矩阵求逆等。

其次，ndarray支持SIMD（Single Instruction Multiple Data），可以进一步提升计算性能。

SIMD 的全称是 Single Instruction Multiple Data，中文名“单指令多数据”。顾名思义，一条指令处理多个数据。

如上图所示

(a) 是普通的标量计算，加法指令一次只能对两个数执行一个加法操作。

(b)是SIMD向量计算，SIMD加法指令一次可以对两个数组（向量）执行加法操作。

ndarray简明教程

看了上面的介绍我猜你已经迫不及待想尝试一下ndarray的威力了。下面我们开始ndarray的学习。

首先让我们引入ndarray包。

传统的Rust程序有Cargo进行包管理，只需要在cargo.toml的[dependencies]中加入

ndarray = "0.15.6"

然后在main.rs开头加入

use ndarray::prelude::*;

即可引入ndarray。prelude涵盖了几乎一切我们需要的模块。

但是在交互环境下没有Cargo，我们需要用evcxr的:dep命令来引入包。在Jupyter中输入如下代码：

:dep ndarray = {version = "0.15.6"}

use ndarray::prelude::*;

创建数组

我们先来看一下如何创建数组：

let arr = array![1., 2., 3., 4., 5., 6.];
println!("1D array: {}", arr);

ndarray提供了array!宏来创建数组。array!宏会自动判断创建那种类型的ArrayBase。上面的示例代码创建了一个1-D类型的ArrayBase，即一维数组。基类ArrayBase实现了std::fmt::Display方法，因此创建的对象可以直接利用println!宏输出。

我们比较一下ndarray与rust自带数组和vec的区别：

let arr1 = array![1., 2., 3., 4., 5., 6.];
println!("1D array: \t{}", arr1);

let ls1 = [1., 2., 3., 4., 5., 6.];
println!("1D list: \t{:?}", ls1);

let vec1 = vec![1., 2., 3., 4., 5., 6.];
println!("1D vector: \t{:?}", vec1);

输出结果如下：

1D array:     [1, 2, 3, 4, 5, 6]
1D list:     [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
1D vector:     [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]

注意，array!将浮点型转成了整形，因为给定的数小数位都为0。

向量加法

下面我们让两个数组中的元素两两相加：

let arr2 = array![1., 2.2, 3.3, 4., 5., 6.];
let arr3 = arr1 + arr2;
println!("1D array: \t{}", arr3);

对比Rust自带数组和vec的实现，你就能发现ndarray多么简单自然。

let arr2 = array![1., 2.2, 3.3, 4., 5., 6.];
let arr3 = arr1 + arr2;
println!("1D array: \t{}", arr3);

let ls2 = [1., 2.2, 3.3, 4., 5., 6.];
let mut ls3 = ls1.clone();
for i in 1..ls2.len(){
    ls3[i] = ls1[i] + ls2[i];
}
println!("1D list: \t{:?}", ls3);

let vec2 = vec![1., 2.2, 3.3, 4., 5., 6.];
let vec3: Vec<f64> = vec1.iter().zip(vec2.iter()).map(|(&e1, &e2)| e1 + e2).collect();
println!("1D vec: \t{:?}", vec3);

输出结果：

1D array:     [2, 4.2, 6.3, 8, 10, 12]
1D list:     [1.0, 4.2, 6.3, 8.0, 10.0, 12.0]
1D vec:     [2.0, 4.2, 6.3, 8.0, 10.0, 12.0]

上面的对比可以发现，Rust自带数组和vec都需要循环或迭代逐个元素相加求值，而ndarray只需要简单的+即可，即简单又直观。对数据和算法工程师来说，ndarray可以极大提高开发效率。

二维数组和高维数组

接下来的演示就不再给大家比较Rust自带数组和vec的实现了。上一节已经充分说明用Rust自带数据类型和数据结构来实现数组操作非常繁琐，而ndarray则十分简洁直观。所以下面我们将聚集ndarray。

ndarray提供多种方法创建和实例化二维数组。看下面代码：

let arr4 = array![[1., 2., 3.], [ 4., 5., 6.]];
let arr5 = Array::from_elem((2, 1), 1.);
let arr6 = arr4 + arr5;
println!("2D array:\n{}", arr6);

输出如下：

2D array:
[[2, 3, 4],
 [5, 6, 7]]

上面代码展示了2种创建二维数组的方法。

• array!宏需要指定数组中的所有元素
• Array::from_elem需要我们给定Shape（在上面案例中是（2，1））和填充元素（在上面案例中是1.0），该方法会用给定元素填充指定大小的二维数组

ndarray还提供了快速创建常用二维数组（矩阵）的方法。看下面代码：

let arr7 =  Array::<f64, _>::zeros(arr6.raw_dim());
let arr8 = arr6 * arr7;
println!("{}", arr8);

输出：

[[0, 0, 0],
 [0, 0, 0]]

Array::zeros(Shape)用0填充指定Shape大小的矩阵。

这里的<f64, _>要解释一下，因为Rust有着严格的类型，编译器有时无法通过0推断出数据类型，所以需要我们显式地告诉编译器这里填充的数据是什么类型。后面的_表示让编译器推断这里的数据类型，此处的数据类型是Shape。

.raw_dim()方法我想你已经猜到了它的作用，没错，就是返回矩阵的维度（dimension）。

在线性代数和机器学习中经常会用到单位矩阵（主对角线元素为1，其余元素为0的二维矩阵）。ndarray提供了专门的eye方法来生成单位矩阵。请看下面的代码：

let identity: Array2<f64> = Array::eye(3);
println!("\n{}", identity);

输出：

[[1, 0, 0],
 [0, 1, 0],
 [0, 0, 1]]

eye方法需要我们传入矩阵大小，因为单位矩阵一定是方阵，所以只需要传一个整数即可。这里要注意的是identity的类型我们明确定义为Array2类型，Array2是ArrayBase子类，表示一个二位数组。

单位矩阵（数学上一般用$I$表示）有很多性质，比如任何矩阵乘以单位矩阵都等于它本身： $A · I = A = I · A$，就像我们经典代数中的1一样。我们可以验证一下单位矩阵的这个性质。

let arr9 = array![[1., 2., 3.], [ 4., 5., 6.], [7., 8., 9.]];
let arr10 = arr9 * identity;
println!("{}", arr10);

结果输出：

[[1, 0, 0],
 [0, 5, 0],
 [0, 0, 9]]

输出结果似乎不是arr9本身。这是因为矩阵运算中*和·不是一回事。*也叫标量乘法，是将两个矩阵对应位置元素相乘。而·是将前一个矩阵的行向量与后一个矩阵的列向量相乘。因此·对矩阵的维度有要求，且不满足交换律。相比*,·乘在机器学习中用的更多，ndarray提供了dot()函数来计算·乘。

let arr11: Array2<f64> = arr9.dot(&identity);
println!("{}", arr11);

这次输出就跟我们的预期一致了：

[[1, 2, 3],
 [4, 5, 6],
 [7, 8, 9]]

同样的道理，我们可以声明3维度或更高维度的数组。例如：

let arr12 = Array::<i8, _>::ones((2, 3, 2, 2));
println!("MULTIDIMENSIONAL\n{}", arr12);

输出的是一个$2 \times 3 \times 2 \times 2$数组

MULTIDIMENSIONAL
[[[[1, 1],
   [1, 1]],

  [[1, 1],
   [1, 1]],

  [[1, 1],
   [1, 1]]],


 [[[1, 1],
   [1, 1]],

  [[1, 1],
   [1, 1]],

  [[1, 1],
   [1, 1]]]]

用ndarray-rand生成随机数

ndarray-rand提供了比rand更强大的随机功能，且能更好得跟ndarray配合。

要想在交互编程环境中引入ndarray-rand，需要加入如下代码：

:dep ndarray-rand={version = "0.14.0"}

use ndarray_rand::{RandomExt, SamplingStrategy};
use ndarray_rand::rand_distr::Uniform;

下面我们生成一个大小为$2 \times 5$的在满足[1, 10]之间正态分布的矩阵。

let arr13 = Array::random((2, 5), Uniform::new(0., 10.));
println!("{:5.2}", arr13);

输出结果为：

[[ 1.64,  8.13,  5.92,  5.88,  0.45],
 [ 7.13,  5.79,  8.21,  4.64,  6.04]]

上面的代码每次运行输出结果都不同，因为每次都是随机生成的，但生成的随机数一定满足[1, 10]之间的正态分布。

我们还可以实现类似Python中random.choice()的功能，随机从数组中选取元素

let arr14 = array![1., 2., 3., 4., 5., 6.];
let arr15 = arr14.sample_axis(Axis(0), 2, SamplingStrategy::WithoutReplacement);
println!("Sampling from:\t{}\nTwo elements:\t{}", arr14, arr15);

输出结果为：

Sampling from:    [1, 2, 3, 4, 5, 6]
Two elements:    [1, 3]

接下来给大家演示另一种随机采样方法，这种方法需要用到rand包。

首先我们引入需要用到的模块

use ndarray_rand::rand as rand;
use rand::seq::IteratorRandom;

这里我们将ndarray-rand包中的rand重命名为rand。

然后我们就可以仿照rand文档中的示例，用ndarray-rand实现相同的功能

let mut rng = rand::thread_rng();
let faces = "😀😎😐😕😠😢";
let arr16 = Array::from_shape_vec((2, 2), faces.chars().choose_multiple(&mut rng, 4)).unwrap();
println!("Sampling from:\t{}", faces);
println!("Elements:\n{}", arr16);

上面的代码，我们首先创建了一个thread_rng实例，然后创建一个字符串保存一组emoji字符，接下来的代码将从这个字符串中随机选取4个字符，组成一个大小为$2 \times 2$的矩阵。这里随机选取用到了IteratorRandom中的choose_multiple方法，随机从字符串中选择4个字符。

这里注意，choose_multiple的采样数量不能超过字符串中的字符数，否则会报panic。

另外，Array::from_shape_vec返回的是Result类型，表示能否创建给定大小的数组。我们上面的代码简单的unwrap()一下获取Result中的结果。

最终输出结果为：

Sampling from:    😀😎😐😕😠😢
Elements:
[[😠, 😢],
 [😐, 😕]]

用ndarray-stats实现统计

接下来我们看一下ndarray-stats，同时介绍一个经常与ndarray-stats配合使用的包noisy_float。

我们将利用ndarray-stats 来验证一下 ndarray-rand生成的正态分布是否真的是正态分布。

首先，我们生成一个大小为$10000 \times 2$的矩阵，矩阵中的元素值满足标准正态分布。

use ndarray_rand::rand_distr::{Uniform, StandardNormal}; 

let arr17 = Array::<f64, _>::random_using((10000,2), StandardNormal, &mut rand::thread_rng());

接下来，我们引入需要用到的统计模块。首先将ndarray-stats和noisy_float引入交互式编程环境：

:dep ndarray-stats = {version="0.5.1"}
:dep noisy_float = {version="0.2.0"}

然后引入如下模块:

use ndarray_stats::HistogramExt;
use ndarray_stats::histogram::{strategies::Sqrt, GridBuilder};
use noisy_float::types::{N64, n64};

接着我们需要将数组中的元素从float转为noisy float。这里我不过多的解释noisy float，你可以简单的将noisy float理解为就是一个float，但是当赋与其非法值时（如NaN）就会报panic。由于noisy float不存在非法值，所以它是安全可排序的，这是ndarray-stats创建直方图等操作必要的条件。

要对ndarray中每一个元素进行操作，我们可以用mapv()。mapv()跟map()类似，传入一个闭包对每一个迭代值机型操作。

let data = arr17.mapv(|e| n64(e));

将数组中的元素转成noisy float后，我们就可以创建统计直方图了。直方图需要以一定的步长切分数据，在ndarray-stats中称之为grid。ndarray-stats提供了GridBuilder来完成数据切分工作。这里我们使用strategies::Sqrt（很多应用都采用此策略，例如我们熟悉的Excel）来帮我们自动推断最佳的切分步长。

let grid = GridBuilder::<Sqrt<N64>>::from_array(&data).unwrap().build();

有了切分步长后我们就可以将其应用到原始数据中得到直方图数据

let histogram = data.histogram(grid);

我们关注的是每个区间内元素的数量，所以我们可以用counts()方法来统计数量

let histogram_matrix = histogram.counts();

这里得到的histogram_matrix 是一个计数矩阵，统计了每个grid内元素的数量。我们的直方图一般是按列统计的，所以我们需要将histogram_matrix按有轴累加起来，这样得到的一维数组就是直方图的统计计数了。

let data = histogram_matrix.sum_axis(Axis(0));
println!("{}", data);

输出结果为：

[2, 0, 1, 0, 5, 5, 1, 4, 9, 2, 4, 12, 8, 13, 17, 11, 20, 32, 32, 30, 32, 43, 47, 66, 65, 89, 110, 114, 110, 131, 147, 170, 173, 194, 209, 206, 255, 259, 264, 280, 292, 255, 288, 294, 335, 300, 306, 324, 287, 323, 275, 280, 284, 250, 244, 247, 192, 230, 178, 191, 162, 147, 150, 128, 110, 99, 83, 93, 83, 68, 47, 50, 30, 32, 34, 21, 24, 13, 14, 15, 9, 12, 4, 4, 7, 2, 3, 1, 0, 1, 1, 1, 0, 0, 3, 0, 0, 0, 0, 2, 0]

虽然输出结果不如图像直观，但也能感觉到确实是服从正态分布的。为了让大家直观看到结果，我们来将数据可视化。数据可视化是数据分析和机器学习中非常重要的一块，后面我会专门一张讲解如何在交互环境下用plotters进行数据可视化。在这里主要让大家可视化地看一下结果。

:dep plotters = {version="0.3.4", default_features = false, features = ["evcxr", "all_series", "all_elements"]}

extern crate plotters;
use plotters::prelude::*;

let figure = evcxr_figure((640, 480), |root| {
    root.fill(&WHITE);
    let histogram_matrix = histogram.counts();
    let data = histogram_matrix.sum_axis(Axis(0));
    let mut chart = ChartBuilder::on(&root)
        .set_label_area_size(LabelAreaPosition::Left, 40)
        .set_label_area_size(LabelAreaPosition::Bottom, 40)
        .caption("随机数分布图", ("sans-serif", 30))
        .build_cartesian_2d((0..data.len()).into_segmented(), 0..300)
        .unwrap();

    chart
        .configure_mesh()
        .disable_x_mesh()
        .bold_line_style(&WHITE.mix(0.3))
        .axis_desc_style(("sans-serif", 15))
        .draw()?;

    chart.draw_series((0..).zip(data.iter()).map(|(x, y)| {
        let x0 = SegmentValue::Exact(x);
        let x1 = SegmentValue::Exact(x + 1);
        let mut bar = Rectangle::new([(x0, 0), (x1, (*y) as i32)], RED.filled());
        bar.set_margin(0, 0, 1, 1);
        bar
    }))?;

    Ok(())
});
figure

上面代码会将data中的数据以柱状图的形式显示。

结论

本章我们学习了如何使用ndarray以及它与自带数组和vec的区别。我们还学习了ndarray生态中随机和统计相关的包，这些包一定程度上提供了scipy的相关功能。

其实ndarray的功能远远不止这些，基本上numpy有的功能ndarray都有对应的实现。这篇文章详细对比了ndarray 和 numpy的函数功能，熟悉numpy的同学可以对比着学ndarray，效率会非常高。后面我也会专门做一篇numpy和ndarray方法的详细对比。