Rust机器学习之ndarray

Rust机器学习之ndarray

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

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

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

在这里插入图片描述

文章目录

什么是ndarray?

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

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

可以说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::*;

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

但是在交互环境下没有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类型,Array2ArrayBase子类,表示一个二位数组。

单位矩阵(数学上一般用 I I I表示)有很多性质,比如任何矩阵乘以单位矩阵都等于它本身: A ⋅ I = A = I ⋅ A A · I = A = I · A 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 × 3 × 2 × 2 2 \times 3 \times 2 \times 2 2×3×2×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 × 5 2 \times 5 2×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 × 2 2 \times 2 2×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 × 2 10000 \times 2 10000×2的矩阵,矩阵中的元素值满足标准正态分布。

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

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

接下来,我们引入需要用到的统计模块。首先将ndarray-statsnoisy_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创建直方图等操作必要的条件。noisy float的具体介绍请参考noisy_float文档

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

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

将数组中的元素转成noisy float后,我们就可以创建统计直方图了。直方图需要以一定的步长切分数据,在ndarray-stats中称之为gridndarray-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都有对应的实现。这篇文章详细对比了ndarraynumpy的函数功能,熟悉numpy的同学可以对比着学ndarray,效率会非常高。后面我也会专门做一篇numpyndarray方法的详细对比。
转自:https://jarod.blog.csdn.net/article/details/127807622

posted @ 2023-06-05 21:46  MasonLee  阅读(558)  评论(0编辑  收藏  举报