在本篇文章中，作者將在 Rust 上移植一個(gè)簡(jiǎn)單的神經(jīng)網(wǎng)絡(luò)實(shí)現(xiàn)。他的目標(biāo)是探索 Rust 中的數(shù)據(jù)科學(xué)工作流在性能和工程學(xué)上的表現(xiàn)。

Python 實(shí)現(xiàn)

第一章描述了一個(gè)非常簡(jiǎn)單的單層神經(jīng)網(wǎng)絡(luò)。這個(gè)神經(jīng)網(wǎng)絡(luò)可以使用基于隨機(jī)梯度下降的機(jī)器學(xué)習(xí)算法，對(duì)來(lái)自于 MNIST 數(shù)據(jù)集的手寫(xiě)數(shù)字進(jìn)行分類。這聽(tīng)起來(lái)挺復(fù)雜，這些東西也確實(shí)在上世紀(jì) 80 年代中期是最先進(jìn)的，但是實(shí)際上，這全部是由一段 150 行的 Python 代碼做出來(lái)的，而且這些代碼有很多人評(píng)論過(guò)。

如果你已經(jīng)知道了這一節(jié)的內(nèi)容(神經(jīng)網(wǎng)絡(luò)基礎(chǔ)知識(shí))，那么我建議你可以跳過(guò)去，當(dāng)然如果想再?gòu)?fù)習(xí)一下神經(jīng)網(wǎng)絡(luò)的基礎(chǔ)知識(shí)也是可以看這一節(jié)的�；蛘卟灰�只關(guān)注代碼，沒(méi)有必要特別細(xì)致的理解代碼為什么以現(xiàn)有的方式運(yùn)行，而應(yīng)該關(guān)注 Python 方法和 Rust 方法的不同。

代碼中的基礎(chǔ)數(shù)據(jù)容器是一個(gè) Network 類，它表示一個(gè)對(duì)層數(shù)和每層網(wǎng)絡(luò)數(shù)量可控制的神經(jīng)網(wǎng)絡(luò)。在 Network 類的內(nèi)部，用 2D NumPy 數(shù)組組成的列表表示類的數(shù)據(jù)。網(wǎng)絡(luò)各層用一個(gè)二維的權(quán)重?cái)?shù)組和一個(gè)一維的偏差數(shù)組來(lái)表示，這些數(shù)組被包含在叫做 biases 和 weights 的 Network 類的屬性中。這些都是二維數(shù)組列表。biases 屬性是列向量，但是會(huì)利用虛擬維度被存儲(chǔ)成二維的數(shù)組。Network 的構(gòu)造函數(shù)如下：

class Network(object):

def __init__(self, sizes):
"""The list ``sizes`` contains the number of neurons in the
respective layers of the network. For example, if the list
was [2, 3, 1] then it would be a three-layer network, with the
first layer containing 2 neurons, the second layer 3 neurons,
and the third layer 1 neuron. The biases and weights for the
network are initialized randomly, using a Gaussian
distribution with mean 0, and variance 1. Note that the first
layer is assumed to be an input layer, and by convention we
won't set any biases for those neurons, since biases are only
ever used in computing the outputs from later layers."""
self.num_layers = len(sizes)
self.sizes = sizes
self.biases = [np.random.randn(y, 1) for y in sizes[1:]]
self.weights = [np.random.randn(y, x)
for x, y in zip(sizes[:-1], sizes[1:])]

在這個(gè)簡(jiǎn)單的實(shí)現(xiàn)中，屬性 biases 和 weights 通過(guò)描述標(biāo)準(zhǔn)正態(tài)分布來(lái)初始化。正態(tài)分布的均值為 0，標(biāo)準(zhǔn)差為 1。我們還可以清楚地看到，biases 是如何被初始化為列向量的。

Network 類公開(kāi)了兩個(gè)可以被用戶直接調(diào)用的方法。第一個(gè)方法是 evaluate 方法，這個(gè)方法可以通過(guò)網(wǎng)絡(luò)嘗試識(shí)別一系列測(cè)試圖片中的數(shù)字，然后基于先驗(yàn)已知的正確結(jié)果，對(duì)識(shí)別的結(jié)果進(jìn)行打分。第二個(gè)方法是 SGD 方法，這個(gè)方法可以通過(guò)遍歷一組圖片來(lái)執(zhí)行隨機(jī)梯度下降的過(guò)程。這個(gè)過(guò)程包括：將整組圖像分解成小的類別，基于各個(gè)小類別圖像來(lái)更新網(wǎng)絡(luò)狀態(tài)，更新用戶指定的學(xué)習(xí)率，eta，以及在用戶隨機(jī)指定數(shù)量的一系列小類別圖像上重新運(yùn)行以上的訓(xùn)練步驟。每組小分類圖像和網(wǎng)絡(luò)的更新的核心算法如下代碼所示：

def update_mini_batch(self, mini_batch, eta):
"""Update the network's weights and biases by applying
gradient descent using backpropagation to a single mini batch.
The ``mini_batch`` is a list of tuples ``(x, y)``, and ``eta``
is the learning rate."""
nabla_b = [np.zeros(b.shape) for b in self.biases]
nabla_w = [np.zeros(w.shape) for w in self.weights]
for x, y in mini_batch:
delta_nabla_b, delta_nabla_w = self.backprop(x, y)
nabla_b = [nb+dnb for nb, dnb in zip(nabla_b, delta_nabla_b)]
nabla_w = [nw+dnw for nw, dnw in zip(nabla_w, delta_nabla_w)]
self.weights = [w-(eta/len(mini_batch))*nw
for w, nw in zip(self.weights, nabla_w)]
self.biases = [b-(eta/len(mini_batch))*nb
for b, nb in zip(self.biases, nabla_b)]

對(duì)于小類別集合中的每個(gè)經(jīng)受訓(xùn)練的圖像，我們通過(guò)反向傳播算法(在 backprop 函數(shù)中實(shí)現(xiàn))積累成本函數(shù)的梯度估計(jì)。當(dāng)程序跑完了小分類圖像集合，會(huì)根據(jù)估計(jì)梯度調(diào)整權(quán)重和偏差。因?yàn)槲覀兿Ｍ�得到小分類集合中所有估�?jì)的平均梯度，所以要更新的數(shù)據(jù)包括分母中的 len(mini_batch) 值。我們還可以通過(guò)調(diào)整學(xué)習(xí)速率和 eta 值來(lái)控制權(quán)重和偏差的更新速度，從而可以全局調(diào)整每個(gè)小分類集合更新的大小。backprop 函數(shù)開(kāi)始于給定輸入圖像的網(wǎng)絡(luò)預(yù)期，然后通過(guò)網(wǎng)絡(luò)向后面運(yùn)行，以在網(wǎng)絡(luò)中通過(guò)層來(lái)傳播這些錯(cuò)誤，最后計(jì)算出神經(jīng)網(wǎng)絡(luò)的成本梯度。這需要大量的數(shù)據(jù)調(diào)整，這也是我在移植到 Rust 時(shí)，花費(fèi)了大量時(shí)間的地方。不過(guò)我認(rèn)為我深入這塊花費(fèi)了過(guò)于長(zhǎng)的時(shí)間。如果你想更深入的了解細(xì)節(jié)，請(qǐng)看這本書(shū)的第 2 章。

Rust 實(shí)現(xiàn)

第一步要搞清楚如何加載數(shù)據(jù)。這一步太繁瑣了，我寫(xiě)了一篇專門(mén)的文章來(lái)介紹。按照順序，我必須弄清楚如何用 Rust 實(shí)現(xiàn) Python 的 Network 類。最終我決定使用結(jié)構(gòu)體(struct)：

use ndarray::Array2;

#[derive(Debug)]
struct Network {
num_layers: usize,
sizes: Vec<usize>,
biases: Vec<Array2<f64>>,
weights: Vec<Array2<f64>>,
}

和 Python 實(shí)現(xiàn)大致一樣，依據(jù)每層神經(jīng)網(wǎng)絡(luò)的數(shù)量初始化結(jié)構(gòu)體。

use rand::distributions::StandardNormal;
use ndarray::{Array, Array2};
use ndarray_rand::RandomExt;

impl Network {
fn new(sizes: &[usize]) -> Network {
let num_layers = sizes.len();
let mut biases: Vec<Array2<f64>> = Vec::new();
let mut weights: Vec<Array2<f64>> = Vec::new();
for i in 1..num_layers {
biases.push(Array::random((sizes[i], 1), StandardNormal));
weights.push(Array::random((sizes[i], sizes[i - 1]), StandardNormal));
}
Network {
num_layers: num_layers,
sizes: sizes.to_owned(),
biases: biases,
weights: weights,
}
}
}

有一個(gè)和 Python 實(shí)現(xiàn)的不同點(diǎn)。我們?cè)?Python 中，使用 numpy.random.randn 來(lái)初始化權(quán)重和偏差，而在 Rust 中，我們使用 ndarray::Array::random 函數(shù)接受一個(gè) rand::distribution::Distributionas 類型的參數(shù)和一個(gè)其他參數(shù)，并允許選擇任意分布，來(lái)完成初始化。在這種情況下，我們使用了 rand::distributions::StandardNormal 做分布。值得注意的是，這使用了在三個(gè)不同包中定義的接口：其中兩個(gè)接口，ndarray 自身和 ndarray-rand 由 ndarray 的作者維護(hù)，剩下一個(gè)由其它的開(kāi)發(fā)者維護(hù)。

這些統(tǒng)一類包的優(yōu)勢(shì)

原則上，一個(gè)好處是，隨機(jī)數(shù)生成不會(huì)在 ndarray 代碼庫(kù)中被單獨(dú)隔離，如果新的隨機(jī)數(shù)分布或者功能被加到了 rand 中，ndarray 和在 Rust 系統(tǒng)中的類都可以同等的使用。另一方面，需要為了各種裝箱操作，在不同文件之間進(jìn)行引用，而不是可以在一個(gè)集中的地方進(jìn)行查看，這樣增加了一些學(xué)習(xí)成本。我也有個(gè)特殊情況，也算是運(yùn)氣不好，rand 發(fā)布了改變公共 api 的新版本時(shí)，我正在開(kāi)發(fā)這個(gè)工程。這導(dǎo)致了，依賴于 0.6 版本的 rand 的 ndarray-rand 和依賴于 0.7 版本的 rand 的我的工程，之間的不一致。

我了解到 cargo 和 rust 的構(gòu)建系統(tǒng)可以很好地處理這類問(wèn)題，但至少在這種情況下，我遇到了一個(gè)令人困惑的錯(cuò)誤信息。這則錯(cuò)誤信息是關(guān)于我的隨機(jī)分布如何不滿足 Distribution(分布)特征的。然而，我的分布式正確的，它的隨機(jī)分布特征滿足了 0.7 版本的 rand，而不是 ndarray-rand 依賴的 rand0.6 版本。但是因?yàn)檠b箱的版本信息不出現(xiàn)在錯(cuò)誤信息當(dāng)中，所以會(huì)讓人非常的困惑。我最后提交了一個(gè) issue 。我發(fā)現(xiàn)對(duì)于 Rust 語(yǔ)言來(lái)說(shuō)，來(lái)自包裝箱的各種不一致接口的讓人困惑的錯(cuò)誤信息，是一個(gè)長(zhǎng)期存在的問(wèn)題。希望在未來(lái)，Rust 可以產(chǎn)生更多有用的錯(cuò)誤信息。

最后，作為一個(gè)新用戶，這種關(guān)注點(diǎn)的分離給我(的理解上)帶來(lái)了許多阻力。在 Python 中，我可以簡(jiǎn)單的做 import numpy 操作便可以完成導(dǎo)入。我認(rèn)為 NumPy 在完全單片化的方向上走的太遠(yuǎn)了。它最初被編寫(xiě)的時(shí)候，用 C 語(yǔ)言擴(kuò)展的 Python 代碼，在打包和發(fā)布上，要比現(xiàn)在困難的多。我認(rèn)為在一個(gè)極端方向上走的太遠(yuǎn)，會(huì)讓一個(gè)語(yǔ)言或者工具的生態(tài)系統(tǒng)變得很難學(xué)習(xí)。

類型和所有權(quán)

下一步我將詳細(xì)介紹 update_mini_batch 的 Rust 版本：

impl Network {
fn update_mini_batch(
&mut self,
training_data: &[MnistImage],
mini_batch_indices: &[usize],
eta: f64,
) {
let mut nabla_b: Vec<Array2<f64>> = zero_vec_like(&self.biases);
let mut nabla_w: Vec<Array2<f64>> = zero_vec_like(&self.weights);
for i in mini_batch_indices {
let (delta_nabla_b, delta_nabla_w) = self.backprop(&training_data[*i]);
for (nb, dnb) in nabla_b.iter_mut().zip(delta_nabla_b.iter()) {
*nb += dnb;
}
for (nw, dnw) in nabla_w.iter_mut().zip(delta_nabla_w.iter()) {
*nw += dnw;
}
}
let nbatch = mini_batch_indices.len() as f64;
for (w, nw) in self.weights.iter_mut().zip(nabla_w.iter()) {
*w -= &nw.mapv(|x| x * eta / nbatch);
}
for (b, nb) in self.biases.iter_mut().zip(nabla_b.iter()) {
*b -= &nb.mapv(|x| x * eta / nbatch);
}
}
}

該函數(shù)使用了我定義的兩個(gè)簡(jiǎn)短的輔助函數(shù)，使得代碼更簡(jiǎn)潔了一些：

fn to_tuple(inp: &[usize]) -> (usize, usize) {
match inp {
[a, b] => (*a, *b),
_ => panic!(),
}
}

fn zero_vec_like(inp: &[Array2<f64>]) -> Vec<Array2<f64>> {
inp.iter()
.map(|x| Array2::zeros(to_tuple(x.shape())))
.collect()
}

和 Python 的實(shí)現(xiàn)版本相比較，調(diào)用 update_mini_batch 的方式有些不同。沒(méi)有直接傳遞對(duì)象列表，反而傳遞了對(duì)完整集合中全套訓(xùn)練數(shù)據(jù)和一份索引的引用。這樣，更容易理解沒(méi)有觸發(fā)器的借用檢查器。

創(chuàng)建 nabla_b 和 nabla_win zero_vec_like 和我們?cè)?Python 中使用的列表解析非常類似。有一個(gè)挫折讓我有些沮喪，因?yàn)槿绻�我嘗試使用 Array2::zeros 創(chuàng)建一個(gè)用 0 填充的數(shù)組，并將它傳遞給一個(gè)特定形狀的 slice 或 Vec，我就會(huì)得到一個(gè) ArrayD 的實(shí)例。為了獲得 Array2 對(duì)象(顯然這是一個(gè)二維數(shù)組，而不是一個(gè)通用的 D 緯數(shù)組)，我需要向 Array::zeros 傳遞一個(gè)元素。然而，由于 ndarray::shape 返回一個(gè)切片(slice)，我需要使用 to_tuple 函數(shù)，將這個(gè)切片轉(zhuǎn)換為一個(gè)元組。這些事情在 Python 中可以被隱藏，但是在 Rust 中，切片(slice)和元組(tuple)之間的不同(造成的影響)變得非常大，就和在這個(gè) API 中的情況一樣。

通過(guò)反向傳播對(duì)估計(jì)的權(quán)重和偏差進(jìn)行更新的代碼具有和 Python 實(shí)現(xiàn)的版本非常相似的結(jié)構(gòu)。我們?cè)谛》诸愔杏?xùn)練每個(gè)示例圖像，并獲得二次成本梯度的估計(jì)值作為偏差和權(quán)重相關(guān)的一個(gè)函數(shù)：

let (delta_nabla_b, delta_nabla_w) = self.backprop(&training_data[*i]);

然后累積這些估計(jì)值：

for (nb, dnb) in nabla_b.iter_mut().zip(delta_nabla_b.iter()) {
*nb += dnb;
}
for (nw, dnw) in nabla_w.iter_mut().zip(delta_nabla_w.iter()) {
*nw += dnw;
}

等到我們完成了小分類的處理，我們就會(huì)根據(jù)學(xué)習(xí)率更新權(quán)重和偏差。

let nbatch = mini_batch_indices.len() as f64;
for (w, nw) in self.weights.iter_mut().zip(nabla_w.iter()) {
*w -= &nw.mapv(|x| x * eta / nbatch);
}
for (b, nb) in self.biases.iter_mut().zip(nabla_b.iter()) {
*b -= &nb.mapv(|x| x * eta / nbatch);
}

這個(gè)例子說(shuō)明了，Rust 對(duì)數(shù)組數(shù)據(jù)在工程上的處理和 Python 相比有區(qū)別的。首先，我們不用浮點(diǎn)數(shù) eta/nbatch 乘以數(shù)組，而是使用 Array::mapv，并定義一個(gè)內(nèi)連閉包，以便在整個(gè)數(shù)組上以矢量化的方式進(jìn)行映射。在 Python 中，由于方法調(diào)用比較慢，所以這些事情不會(huì)處理的太快。而在 Rust 中，則不會(huì)出現(xiàn)這種情況。當(dāng)我們減去時(shí)，還需要借用帶 & 符號(hào) mapv 的返回值，以免我們?cè)诘�代它時(shí)消耗數(shù)組數(shù)據(jù)。在 Rust 中，需要仔細(xì)考慮函數(shù)是否會(huì)消耗數(shù)據(jù)或者引入引用，這導(dǎo)致在概念上，用 Rust 編寫(xiě)這種代碼，比在 Python 中要求更多。另一方面，我對(duì)我的代碼的正確性并且能夠編譯通過(guò)，有了更高的信心。我不確定的是，我寫(xiě)這段代碼很費(fèi)力的原因，是因?yàn)?Rust 真的更難寫(xiě)，還是因?yàn)槲以?Python 和 Rust 上經(jīng)驗(yàn)的不同。

用 Rust 重寫(xiě)這些代碼，然后一切都會(huì)好起來(lái)

在這里，我留下了一些東西，比我開(kāi)始使用的未經(jīng)優(yōu)化的 Python 版本代碼更快。然而，相比 10 倍或是更快的速度，人們可能更期望從像 Python 這樣的動(dòng)態(tài)解釋性語(yǔ)言轉(zhuǎn)變?yōu)橄?Rust 這樣的編譯性能導(dǎo)向語(yǔ)言，并且我也只觀察到了 2 倍的加速�？梢匀ダ斫庖幌挛覟槭裁匆獪y(cè)量 Rust 語(yǔ)言的性能表現(xiàn)。幸運(yùn)的是，這里有一個(gè)非常方便的項(xiàng)目，可以為 Rust 工程生成火焰圖： flamegraph 。這里添加了一個(gè) flamegraph 的子命令 cargo，因此只需要在包中執(zhí)行 cargo flamegraph 即可運(yùn)行代碼，就會(huì)編寫(xiě)出一個(gè)可以在瀏覽器中執(zhí)行的火焰圖 svg 文件(原圖為可交互的 svg 腳本，如果希望嘗試，可以查看原網(wǎng)頁(yè))。

 

 

如果你之前還沒(méi)有看過(guò)一個(gè)火焰圖，(我解釋下)，在例程中發(fā)生的程序運(yùn)行時(shí)間與與該例程的條形寬度成正比。主函數(shù)位于圖形的底部，主函數(shù)調(diào)用的函數(shù)在圖形的頂部。這樣你就可以簡(jiǎn)單查看哪些函數(shù)占用了程序中最多的時(shí)間。圖中非常寬的東西代表了花費(fèi)最多時(shí)間的地方。在調(diào)用棧中非常高和寬的函數(shù)，在代碼上花費(fèi)了大量的時(shí)間�？匆幌律厦娴幕鹧鎴D，我們可以發(fā)現(xiàn)一般的時(shí)間都花費(fèi)在了像名字叫 dgemm_kernel_HASWELL 的這類函數(shù)身上，這類函數(shù)是 OpenBLAS 的線性代數(shù)類庫(kù)，剩下的時(shí)間，花費(fèi)在 update_mini_batch 的數(shù)組和分配數(shù)組之間的添加上。我程序的其它所有部分，對(duì)運(yùn)行時(shí)間的貢獻(xiàn)可以忽略不計(jì)。

如果我們?yōu)?Python 代碼做一個(gè)類似的火焰圖，我們會(huì)發(fā)現(xiàn)一個(gè)相似的情況：大部分時(shí)間花費(fèi)在了做線性代數(shù)函數(shù)上去(在 np.dot 反向傳播例程中調(diào)用的地方)。因此，由于不管是 Rust 還是 Python 花費(fèi)的時(shí)間大部分都在數(shù)值性的線性代數(shù)庫(kù)上，我們就不能夠希望得到一個(gè) 10 倍加速的結(jié)果。

實(shí)際情況比這更糟糕。這本書(shū)中的一個(gè)練習(xí)是重寫(xiě)了使用向量化矩陣乘法的 Python 代碼。在這個(gè)方法中，每個(gè)小分類中的所有樣例的反向傳播發(fā)生在單組矢量化矩陣乘法運(yùn)算中。這需要能夠在 3 維和 2 維數(shù)組之間進(jìn)行矩陣乘法。由于每個(gè)矩陣乘法運(yùn)算使用的數(shù)據(jù)量大于非向量化的情況，OpenBLAS 可以更有效地使用 CPU 緩存和寄存器，基本上可以更好地利用我筆記本電腦上的可用 CPU 資源。重寫(xiě)的 Python 版本要比 Rust 版本更快，又快了大約兩倍左右。

理論上，可以對(duì) Rust 代碼進(jìn)行相同的優(yōu)化。但是對(duì)于高于 2 維(的矩陣)的情況，ndarraycrate還不支持矩陣乘法。也可以使用像 rayon 這樣的庫(kù)在小批量更新上使用線程并行化。我在我的筆記本上嘗試這個(gè)(并行化)沒(méi)有看到任何的加速，但是可能在具有更多 CPU 線程的更強(qiáng)大的機(jī)器上會(huì)有作用。我也可以嘗試使用一個(gè)不同的線性代數(shù)函數(shù)實(shí)現(xiàn)，例如，有 TensorFlow 和 Torch 的 Rust 構(gòu)建，但是在這種情況下，我覺(jué)得我也可以使用那些庫(kù)的 Python 構(gòu)建。

Rust 是否適合數(shù)據(jù)科學(xué)的工作流?

現(xiàn)在我不得不說(shuō)，答案是”未知“。在未來(lái)，當(dāng)我需要編寫(xiě)具有小依賴性的低級(jí)別優(yōu)化代碼時(shí)，我肯定會(huì)使用 Rust。但是，如果把 Rust 作為 Python 和 C++ 的完全替代品，還需要一個(gè)更穩(wěn)定和完善的類庫(kù)生態(tài)系統(tǒng)。

標(biāo)簽： 數(shù)據(jù) 蒲Я煊

版權(quán)申明：本站文章部分自網(wǎng)絡(luò)，如有侵權(quán)，請(qǐng)聯(lián)系：west999com@outlook.com
特別注意：本站所有轉(zhuǎn)載文章言論不代表本站觀點(diǎn)！
本站所提供的圖片等素材，版權(quán)歸原作者所有，如需使用，請(qǐng)與原作者聯(lián)系。