Go科学计算之向量篇 - Gonum实现

RjlPKM

虽无用，但有功。

出于某种目的，有时候我们会将简单的事物复杂化，但其复杂程度不应该超出多数人的理解范围。本文原本试图分别对 Gonum 和 Gosl 两个主要 Go 语言科学计算包关于向量的实现及其基本运算进行整理分析。但写到最后发现篇幅过长，因此将其拆分为 Gonum 篇和 Gosl 篇两部分内容，本文为 Gonum 篇。

‌申明 ：

由于时间和精力有限，本文缺失文献整理部分，如果有侵权之嫌，望及时告知；
由于能力有限，请务必对本文持批判态度，尽信书不如无书。

预备知识

为了便于理解，在正式探究Gonum中的向量类型之前我们先简单的回顾一下 slice 和 struct 两种内置的类型。通过将两者进行简单的组合，可以方便的实现几乎所有我们需要用到的特殊类型。

切片类型

切片类型可以简单的视作简化版的动态数组，而在 Go 中是一个由固定长度的特定类型元素组成的序列，一个数组可以由零个或多个元素组成。我们先看看切片的结构定义， reflect.SliceHeader :

type SlinceHeader struct {
    Data uintptr
    Len int
    Cap int
..

其中，内置的 len 函数返回切片中有效元素的长度，内置的 ‌cap 函数返回切片容量大小，容量必须大于或等于切片的长度。切片的申明主要有以下两种方法：

o := []float64{1, 2, 3, 4} // 方法1
o := make([]float64, 4)    // 方法2

其中方法2中返回值为切片 o 的引用，通过函数传递参数之后可以直接修改原始值。此外可以采用内置函数 new 生成切片，不过一般不推荐。关于切片的进一步说明及相关操作可以在网上查找相应的资料。

结构体类型

Go 通过把使用各种数据类型定义的不同变量组合起来的高级数据类型定义为结构体。一个带有属性的结构体试图表示一个现实世界中的实例，可以采用以下方法简单的定义一个结构体

type identifier struct {
    field1 type1
    field2 type2
}

我们可以采用 type T struct {a, b int} 定义一个简单的结构体。

由于数值计算通常需要存储大量的数据，采用指针进行参数传递无疑可以显著的提高计算效率。我们可以使用 new 函数给一个新的结构体变量分配内存，它返回指向已分配内存的指针： var t *T = new(T) 。进一步，可以使用工厂方法来私有化类型，该方法同时适用于切片和结构体。

package matrix
type matrix struct {
func NewMatrix(params) *matrix {
    o := new(matrix) // 初始化 o
    return o
}

在其他包里使用工厂方法：

package main
import "matrix"
o := matrix.NewMatrix(...)

通过绑定方法，结构体可以作为类的一种简单实现。定义方法的一般格式如下：

func (recv receiver_type) methodName(parameter_list) (return_value_list) { ... }

注：这部分内容主要参考《the-way-to-go_ZH_CN》，建议阅读原文。

向量类型的实现及其基本运算

从实用的角度来说，我们可以直接采用切片来表示向量，并且实现其基本运算。然而，向量作为科学计算包中的核心类型之一，其实现需要同时满足高效且易于扩展的要求，因此其定义要相对复杂一些。在正式介绍 Gonum 和 Gosl 对于向量的实现之前，我打算先简单的介绍一下 python 中 numpy 包对于向量的实现。我将借此来表达对 Google 的强烈不满，在Web领域玩的风生水起的同时请人文关怀一下我们工科生！

Numpy 包采用一维数组简单的表示向量，并且提供了大量的功能函数。向量的一般定义如下：

import numpy as np
a = np.array([0, 1, 2, 3])
b = np.array([1, 2, 3, 4])
c = a + b // 加法运算
d = a - b // 减法运算
e = a * b // 乘法运算 （dot product）
...

写到这里我突然不想写了，我开始无比的怀念 python 对我的溺爱。不过为了证明Go也可以成为数值计算中的一门主力编程语言，我只能抹干眼泪继续前行。

Gonom 主要面向科学计算，试图发展成为类似于 Numpy 和 Scipy 一样完备的科学计算库。为了简化写作，本文仅针对实数向量。

Gonum 封装了 BLAS（Basic Linear Algebra Subprograms，基础线性代数程序集）等久经考验的科学计算库来实现其功能，其中 BLAS 是一个应用程序接口（API）标准，用以规范发布基础线性代数操作的数值库（如矢量或矩阵乘法）。 Gonum 将向量简单的定义为：

// VecDense represents a column vector.
type VecDense struct {
    mat blas64.Vector
    // A BLAS vector can have a negative increment, but allowing this
    // in the mat type complicates a lot of code, and doesn't gain anything.
    // VecDense must have positive increment in this package.
}

其中，Dense 用于区分系数矩阵。 VecDense 的工厂方法实现为：

// NewVecDense creates a new VecDense of length n. If data == nil,
// a new slice is allocated for the backing slice. If len(data) == n, data is
// used as the backing slice, and changes to the elements of the returned VecDense
// will be reflected in data. If neither of these is true, NewVecDense will panic.
// NewVecDense will panic if n is zero.
func NewVecDense(n int, data []float64) *VecDense {
    if n <= 0 {
        if n == 0 {
            panic(ErrZeroLength)
        panic("mat: negative dimension")
    if len(data) != n && data != nil {
        panic(ErrShape)
    if data == nil {
        data = make([]float64, n)
    return &VecDense{
        mat: blas64.Vector{
            N:    n,
            Inc:  1,
            Data: data,
}

其中， n 为向量长度， data 为值， ‌Inc 为步长。值得注意的是，类似于 Numpy， Gunum 并不明显区分向量和矩阵。

在其他包里实用工厂方法：

package main
import (
    "fmt"
    "gonum.org/v1/gonum/mat"
func main() {
    // Initialize with the length of the vector, 
    // followed by a slice of floats containing the data.
    u := mat.NewVecDense(3, []float64{1, 2, 3})
    v := mat.NewVecDense(3, []float64{4, 5, 6})
    fmt.Println("u: ", u)
    fmt.Println("v: ", v)
// output
u:  &{{3 [1 2 3] 1}}
v:  &{{3 [4 5 6] 1}}

我们可以采用 o.len() 来方便的计算向量的长度，(o 代表 VecDense 对象）。向量特定位置的值不能简单采用索引的方式来读取，即 o[i] ，Gonum 分别提供了 o.AtVec 和 o.At 两种方法来实现该目的。其中 o.At 同样适用于矩阵，其定义为：

func (v *VecDense) At(i, j int) float64

对于向量， j 取 0，即通过 o.At(i, 0) 来读取特定位置的值。方法 o.AtVec 仅适用于向量，其定义为：

func (v *VecDense) AtVec(i int) float64

同样，我们需要采用 o.SetVec 方法来修改向量特定位置的值，其定义为：

func (v *VecDense) SetVec(i int, val float64)

在对向量类型进行简单的介绍之后，我们主要介绍关于向量的一些基本运算：

‌加法 ：向量加法可以用方法 o.AddVec 加以实现，其定义为：

func (v *VecDense) AddVec(a, b Vector)

根据该定义我们可以看到该方法并不返回一个新的对象，为了更好的了解其作用机理，我们给出其具体实现如下：

// AddVec adds the vectors a and b, placing the result in the receiver.
func (v *VecDense) AddVec(a, b Vector) {
    ar := a.Len()
    br := b.Len()
    if ar != br {
        panic(ErrShape)
    v.reuseAs(ar)
    aU, _ := untranspose(a)
    bU, _ := untranspose(b)
    if arv, ok := aU.(RawVectorer); ok {
        if brv, ok := bU.(RawVectorer); ok {
            amat := arv.RawVector()
            bmat := brv.RawVector()
            if v != a {
                v.checkOverlap(amat)
            if v != b {
                v.checkOverlap(bmat)
            if v.mat.Inc == 1 && amat.Inc == 1 && bmat.Inc == 1 {
                // Fast path for a common case.
                f64.AxpyUnitaryTo(v.mat.Data, 1, bmat.Data, amat.Data)
                return
            f64.AxpyIncTo(v.mat.Data, uintptr(v.mat.Inc), 0,
                1, bmat.Data, amat.Data,
                uintptr(ar), uintptr(bmat.Inc), uintptr(amat.Inc), 0, 0)
            return
    for i := 0; i < ar; i++ {
        v.setVec(i, a.AtVec(i)+b.AtVec(i))
}

该实现表明为了进行向量的加法运算，我们需要首先初始化一个对象来存放最终的计算结果（感觉有点绕），具体做法如下：

w := mat.NewVecDense(3, nil)
w.AddVec(u, v)

当然，我们也可以选择对其中的一个向量的值用计算得到的结果进行覆盖。

‌减法 ：向量的减法运算实现方法为 o.SubVec ，其使用方法于 o.AddVec 一致，可以定义为：

func (v *VecDense) SubVec(a, b Vector)

‌标量乘法 向量与常数的乘积可以采用 o.ScaleVec 方法，其定义为：

func (v *VecDense) ScaleVec(alpha float64, a Vector)

‌元素积 两个向量的元素积（element-wise product）可以采用 o.MulElemVec 方法计算，其定义为：

func (v *VecDense) MulElemVec(a, b Vector)

点积：值得注意的是，Gonum/mat 模块提供了一个简单的函数来实现向量的点积（dot product），即 mat.Dot ，其定义为：

func Dot(a, b Vector) float64

另外 Gonum/mat 内置了矩阵与向量的乘法，该运算在有限单元法中计算雅可比矩阵时十分有用，其定义为：

func (v *VecDense) MulVec(a Matrix, b Vector)

需要提醒的是，在 Numpy 中需要格外小心该操作，要注意得到的结果是列向量还是行向量。我曾在这上面吃过大亏，感兴趣的请自己手动尝试。

为了直观的认识上述方法，我们给出如下算例：

package main
import (
    "fmt"
    "gonum.org/v1/gonum/mat"
func main() {
    u := mat.NewVecDense(3, []float64{1, 2, 3})
    fmt.Println("u: ", u)
    v := mat.NewVecDense(3, []float64{4, 5, 6})
    fmt.Println("v: ", v)
    w := mat.NewVecDense(3, nil)
    w.AddVec(u, v)
    fmt.Println("u + v: ", w)
    // Add u + alpha * v for some scalar alpha
    w.AddScaledVec(u, 2, v)
    fmt.Println("u + 2 * v: ", w)
    // Subtract v from u
    w.SubVec(u, v)
    fmt.Println("v - u: ", w)
    // Scale u by alpha
    w.ScaleVec(23, u)
    fmt.Println("u * 23: ", w)
    // Compute the dot product of u and v
    // Since float64’s don’t have a dot method, this is not done
    //inplace
    d := mat.Dot(u, v)
    fmt.Println("u dot v: ", d)
    // element-wise product
    w.MulElemVec(u, v)
    fmt.Println("u element-wise product v: ", w)
    // Find length of v
    l := v.Len()
    fmt.Println("Length of v: ", l)
// output
u:  &{{3 [1 2 3] 1}}
v:  &{{3 [4 5 6] 1}}
u + v:  &{{3 [5 7 9] 1}}