纯 vs. 非纯

Ndarray 模块包含许多用于操作和执行多维数组的数学运算的函数。 纯函数（也称为不可变函数）是指那些不修改传入变量但始终返回新变量作为结果的函数。相比之下，非纯函数（也称为可变函数）是指那些在原地修改传入变量的函数。

关于纯函数和非纯函数之间的争论永远不会结束。通常，函数式编程倡导使用不可变数据结构。使用非纯函数使得难以推理代码的正确性，因此在决定使用它们时需要仔细考虑。另一方面，仅因为修改了一个元素就每次生成一个新的 1000 x 1000 矩阵似乎也不太实际。

Owl 引入非纯函数是经过许多慎重和实际的考虑的。使用原地修改的主要动机之一是避免昂贵的内存分配和释放操作，这可以在大型 ndarray 和矩阵涉及时显着提高数值应用程序的运行时性能。

我们是否可以拥有两全其美的最佳部分，即编写功能代码并同时在内存上高效？随着您对 Owl 的了解越来越多，您将意识到通过使用计算图懒惰地评估数学表达式，这是可以实现的。程序员关注功能代码，Owl 的计算图模块负责“危险任务”——高效地分配和管理内存。

Ndarray 模块中的许多纯函数都有其相应的非纯版本，区别在于非纯版本的函数名称末尾有一个额外的下划线 “_”。例如，以下函数是 Arr 模块中的纯函数。

  Arr.sin;;
  Arr.cos;;
  Arr.log;;
  Arr.abs;;
  Arr.add;;
  Arr.mul;;

它们对应的非纯函数如下。

  Arr.sin_;;
  Arr.cos_;;
  Arr.log_;;
  Arr.abs_;;
  Arr.add_;;
  Arr.mul_;;

对于Arr.sin x等一元运算符，情况相当直接，x 将在原地修改。但是，对于Arr.add_scalar_ x a和Arr.add_ x y等二元运算符，情况需要一些澄清。对于Arr.add_scalar x a，x 将在原地修改并存储最终结果，这是微不足道的，因为a是标量。

对于Arr.add_ x y，问题是在两个输入都是 ndarray 时存储最终结果的位置。让我们看一下Arr.add_函数的类型。

val Arr.add_ : ?out:Arr.arr -> Arr.arr -> Arr.arr -> unit

从函数类型中可以看出，输出可以由可选的 out 参数指定。如果输入中缺少 out，那么 Owl 将尝试使用第一个操作数（即 x）来存储最终结果。因为 Owl 中的二元操作符默认支持广播操作，这进一步说明了在使用非纯函数时，第一个参数 x 的每个维度都不能小于第二个参数 y 的相应维度。换句话说，非纯函数只允许将较小的 y 广播到足够大以容纳结果的 x。

Owl 中的大多数二元数学函数都与一个简写操作符相关，例如+，-，*和/。非纯版本也有它们自己的操作符。例如，与Arr.(x + y)对应，它返回一个新的 ndarray 中的结果，您可以编写Arr.(x += y)，该操作将x和y相加并将结果保存到x中。

Ndarray vs. Scalar

有三种类型的 ndarray 操作：map、scan 和 reduce。许多函数可以归类为减少操作，例如 Arr.sum、Arr.prod、Arr.min、Arr.mean、Arr.std 等。Owl 中的所有减少函数都有一个名为 axis 的参数。当在多维数组上应用这些减少操作时，有两种可能的情况：

• 如果明确指定了轴，则 Owl 沿指定的轴减少；
• 如果未指定轴，则 Owl 首先将 ndarray 扁平化为矢量，然后沿轴 0 减少所有元素。

如果传入的 ndarray 已经是一维的，则这两种情况是等效的。在以下代码片段中，a 的形状为 [|3;1;3|]，而 b 的形状为 [|1|]，因为它只包含一个元素。

    let x = Arr.sequential [|3;3;3|];;
    let a = Arr.sum ~axis:1 x;;
    let b = Arr.sum x;;
  

如果您想将 b 中的结果与另一个浮点数相加，您需要通过调用 get 函数检索该值。

    let c = Arr.get b [|0|] in
    c +. 10.;;
  

如果我们总是需要从减少操作的返回值中提取标量值，这看起来并不太方便。对于像 Python 和 Julia 这样的语言来说，这不是问题，因为返回类型是动态确定的。但是，对于 OCaml，这变得具有挑战性：我们要么使用统一的类型，要么实现另一组函数。最终，我们在 Owl 的设计中选择了后者。每个减少操作都有两个版本：

• 一个允许您沿指定轴减少或减少所有元素，但始终返回一个 ndarray；
• 一个仅减少所有元素并始终返回标量值。

两者之间的区别在于返回标量的函数的名称末尾带有额外的撇号 “'” 字符。例如，对于返回 ndarray 的第一种函数类型，它们的函数名称看起来像这样。

    Arr.sum;;
    Arr.min;;
    Arr.prod;;
    Arr.mean;;
    Arr.std;;
  

对于返回标量的第二种函数类型，它们的名称看起来像这样。

    Arr.sum';;
    Arr.min';;
    Arr.prod';;
    Arr.mean';;
    Arr.std';;
  

从技术上讲，Arr.sum' 等效于以下代码。

    let sum' x =
      let y = Arr.sum x in
      Arr.get y [|0|]
  

让我们扩展上一个代码片段，并在 OCaml 的 toplevel 中测试它。然后您将立即了解到两者之间的区别。

    let x = Arr.sequential [|3;3;3|];;
    let a = Arr.sum ~axis:1 x;;
    let b = Arr.sum x;;
    let c = Arr.sum' x;;
  

规则和约定通常揭示了设计中的权衡。通过澄清限制，我们希望程序员能够在特定情境中选择正确的函数使用。

Infix Operators

Owl 中的运算符是在 Owl_operator 模块中定义的函数子中实现的。这些运算符分为 Basic、Extend、Matrix 和 Ndarray 四个模块类型签名，因为某些操作仅对特定的数据结构有意义。例如，矩阵乘法 *@ 仅在 Matrix 签名中定义。

只要模块实现了模块签名中定义的所有函数，您就可以使用这些函数子生成相应的运算符。在大多数情况下，您不需要直接使用这些函数子在 Owl 中工作，因为我已经为您完成了生成部分。

这些运算符已经包含在每个 Ndarray 和 Matrix 模块中。以下表总结了当前实现的运算符。在表中，x 和 y 都表示矩阵或 ndarray，而 a 表示标量值。

以下是值得注意的一些事项。

• * 用于逐元素相乘；*@ 用于矩阵相乘。如果阅读 Algodiff 模块的源代码，您可以轻松理解为什么使用 * 进行逐元素相乘（对于矩阵）导致了一致的算法导数实现。

• +$ 有其相应的运算符 $+，如果我们翻转参数的顺序。但是，请非常注意运算符的优先级，因为 OCaml 根据中缀的第一个字符确定优先级。 +$ 保留优先级，而 $+ 不保留。因此，我建议非常谨慎地使用 $+。请始终使用括号显式指定优先级。对于 $-、$* 和 $/ 也是一样。

• 对于比较运算符，例如 = 和 =. 都比较两个变量 x 和 y 中的所有元素。区别在于 = 返回一个布尔值，而 =. 返回形状和类型与 x 和 y 相同的矩阵或 ndarray。在返回的结果中，如果 x 和 y 中相应位置的值满足谓词，则该位置的值为 1，否则为 0。

• 对于以 $ 结尾的比较运算符，它们用于将矩阵/ndarray 与标量值进行比较。

运算符易于使用，以下是一些示例。

let x = Mat.uniform 5 5;;
let y = Mat.uniform 5 5;;

Mat.(x + y);;
Mat.(x * y);;
Mat.(x ** y);;
Mat.(x *@ y);;

(* 比较以下两个返回的结果 *)

Mat.(x > y);;
Mat.(x >. y);;

以下是第一个示例的返回。

Mat.(x > y);;
>- : bool = false

以下是第二个示例的返回。

Mat.(x >. y);;
>- : (float, float64_elt) Owl_dense_matrix_generic.t =
>
>   C0 C1 C2 C3 C4
>R0  0  1  1  0  1
>R1  0  1  1  0  1
>R2  0  0  0  0  0
>R3  1  0  0  0  1
>R4  0  1  0  0  1

现在我相信您已经了解了 > 和 >. 之间的区别，其他二进制比较运算符也是一样的。

请注意，上表未包含扩展的索引和切片运算符，但您可以在 索引和切片章节 中找到详细说明。

运算符扩展

正如您所看到的，上面的运算符不允许在不同的数字类型之间进行互操作（实际上这可能并不是坏事）。例如，除非在 Generic 模块中显式调用 cast 函数，否则您不能将 float32 矩阵添加到 float64 矩阵。对于喜欢 Python 风格的人，Owl.Ext 模块专门为此目的而设计，以使原型制作更快、更容易。一旦打开模块，Ext 就会立即为您提供一组运算符，以便您可以在不同的数字类型上进行互操作，如下所示。如果需要，它会自动为您转换类型。

您可能已经注意到，以 $ 结尾的运算符（例如 +$、-$ 等）在表中消失了，这仅仅是因为我们可以直接将标量与矩阵相加/相减/相乘/相除，因此不再需要这些运算符。对于比较运算符也是类似的，因为我们可以使用相同的 > 运算符来比较矩阵与另一个矩阵，或者将矩阵与标量进行比较，因此我们不再需要 >$。允许互操作使运算符表格变得更短。

目前，Ext 模块中的运算符仅支持密集结构上的互操作。除了二进制运算符之外，Ext 还实现了大多数常见的数学函数，这些函数可以应用于浮点数、复数、矩阵和 ndarray。这些函数包括：

im、re、conj、abs、abs2、neg、reci、signum、sqr、sqrt、cbrt、exp、exp2、expm1、log、log10、log2、log1p、sin、cos、tan、asin、acos、atan、sinh、cosh、tanh、asinh、acosh、atanh、floor、ceil、round、trunc、erf、erfc、logistic、relu、softplus、softsign、softmax、sigmoid、log_sum_exp、l1norm、l2norm、l2norm_sqr、inv、trace、sum、prod、min、max、minmax、min_i、max_i、minmax_i。

请注意，Ext 包含自己的 Ext.Dense 模块，其中进一步包含以下子模块。

• Ext.Dense.Ndarray.S
• Ext.Dense.Ndarray.D
• Ext.Dense.Ndarray.C
• Ext.Dense.Ndarray.Z
• Ext.Dense.Matrix.S
• Ext.Dense.Matrix.D
• Ext.Dense.Matrix.C
• Ext.Dense.Matrix.Z

这些模块只是 Owl.Dense 模块中原始模块的包装器，因此它们已经提供了大多数已实现的 API。这些包装器模块额外提供的是自动为您打包和解包原始数字类型。但是，您当然可以使用原始数据类型，然后使用 Owl_ext_types 中定义的构造函数自己进行包装。构造函数定义如下。

  type ext_typ =
    F   of float
    C   of Complex.t
    DMS of dms
    DMD of dmd
    DMC of dmc
    DMZ of dmz
    DAS of das
    DAD of dad
    DAC of dac
    DAZ of daz
    SMS of sms
    SMD of smd
    SMC of sms
    SMZ of smd
    SAS of sas
    SAD of sad
    SAC of sac
    SAZ of saz

还有相应的 packing 和 unpacking 函数可供使用，请阅读 owl_ext_types.ml <https://github.com/owlbarn/owl/blob/master/src/owl/ext/owl_ext_types.ml>_ 了解更多细节。

让我们通过一些示例了解使用 Ext 模块有多方便。

open Owl.Ext;;

let x = Dense.Matrix.S.uniform 5 5;;
let y = Dense.Matrix.C.uniform 5 5;;
let z = Dense.Matrix.D.uniform 5 5;;

x + F 5.;;
x * C Complex.({re = 2.; im = 3.});;
x - y;;
x / y;;
x *@ y;;

(** ... *)

x > z;;
x >. z;;
(x >. z) * x;;
(x >. F 0.5) * x;;
(F 10. * x) + y *@ z;;

(** ... *)

round (F 10. * (x *@ z));;
sin (F 5.) * cos (x + z);;
tanh (x * F 10. - z);;

(** ... *)

在完成本章之前，我想指出一个注意事项。 Ext 试图通过统一类型来模仿像 Python 这样的动态语言。这阻止了 OCaml 编译器在编译阶段进行类型检查，并在调用函数时引入了额外的开销。因此，除了在顶层进行快速实验之外，我不建议在生产代码中使用 Ext 模块。

模块结构

在 Owl 中，Dense 模块包含密集数据结构的模块。例如，Dense.Matrix 支持密集矩阵的操作。类似地，Sparse 模块包含稀疏数据结构的模块。

Dense.Ndarray;;   (* 密集 ndarray *)
Dense.Matrix;;    (* 密集矩阵 *)

Sparse.Ndarray;;  (* 稀疏 ndarray *)
Sparse.Matrix;;   (* 稀疏 ndarray *)

所有这四个模块都由五个子模块组成，以处理不同类型的数字。

• S 模块支持单精度浮点数 float32；
• D 模块支持双精度浮点数 float64；
• C 模块支持单精度复数 complex32；
• Z 模块支持双精度复数 complex64；
• Generic 模块通过 GADT 支持所有上述数字类型。

使用 Dense.Ndarray，您可以创建一个不超过 16 维的密集 n 维数组。这个限制来源于底层的 Bigarray.Genarray 模块。在实践中，随着维数的增加，空间需求会呈指数级增长，因此这个约束是有意义的。如果您需要的维数超过 16 维，您需要使用 Sparse.Ndarray 创建稀疏数据结构。

数字和精度

在确定了合适的数据结构（密集或稀疏）之后，您可以使用模块中的创建函数（例如 empty、create、zeros、ones 等）创建 ndarray/matrix。数字的类型（实数或复数）及其精度（单精度或双精度）需要作为参数传递给创建函数。

Dense.Ndarray.Generic.zeros Float64 [|5;5|];;
>- : (float, float64_elt) Dense.Ndarray.Generic.t =
>
>   C0 C1 C2 C3 C4
>R0  0  0  0  0  0
>R1  0  0  0  0  0
>R2  0  0  0  0  0
>R3  0  0  0  0  0
>R4  0  0  0  0  0

使用 zeros 函数，创建的数据结构中的所有元素将被初始化为零。

从技术上讲，S、D、C 和 Z 都是带有显式类型信息的 Generic 模块的包装器。因此，如果直接使用这些子模块，则可以省略传递给 Generic 模块的类型构造函数。

Dense.Ndarray.S.zeros [|5;5|];;    (* 单精度实数 ndarray *)
Dense.Ndarray.D.zeros [|5;5|];;    (* 双精度实数 ndarray *)
Dense.Ndarray.C.zeros [|5;5|];;    (* 单精度复数 ndarray *)
Dense.Ndarray.Z.zeros [|5;5|];;    (* 双精度复数 ndarray *)

以下示例是用于密集矩阵的。

Dense.Matrix.S.zeros 5 5;;     (* 单精度实数矩阵 *)
Dense.Matrix.D.zeros 5 5;;     (* 双精度实数矩阵 *)
Dense.Matrix.C.zeros 5 5;;     (* 单精度复数矩阵 *)
Dense.Matrix.Z.zeros 5 5;;     (* 双精度复数矩阵 *)

以下示例是用于稀疏 ndarray 的。

Sparse.Ndarray.S.zeros [|5;5|];;    (* 单精度实数 ndarray *)
Sparse.Ndarray.D.zeros [|5;5|];;    (* 双精度实数 ndarray *)
Sparse.Ndarray.C.zeros [|5;5|];;    (* 单精度复数 ndarray *)
Sparse.Ndarray.Z.zeros [|5;5|];;    (* 双精度复数 ndarray *)

以下示例是用于稀 疏矩阵的。

Sparse.Matrix.S.zeros 5 5;;     (* 单精度实数矩阵 *)
Sparse.Matrix.D.zeros 5 5;;     (* 双精度实数矩阵 *)
Sparse.Matrix.C.zeros 5 5;;     (* 单精度复数矩阵 *)
Sparse.Matrix.Z.zeros 5 5;;     (* 双精度复数矩阵 *)

简而言之，Generic 模块可以执行子模块可以执行的所有操作，但对于某些函数（例如创建函数），您需要显式传递类型信息。

多态函数

在 Generic 模块中，通过模式匹配和 GADT 实现了多态性。这意味着 Generic 模块中的许多函数可以处理前述四种不同类型的数字。

在下面的示例中，我以 Dense.Matrix.Generic 模块中的 sum 函数为例。 sum 函数返回矩阵中所有元素的总和。

open Owl;;

let x = Dense.Matrix.S.eye 5 in
  Dense.Matrix.Generic.sum x;;

let x = Dense.Matrix.D.eye 5 in
  Dense.Matrix.Generic.sum x;;

let x = Dense.Matrix.C.eye 5 in
  Dense.Matrix.Generic.sum x;;

let x = Dense.Matrix.Z.eye 5 in
  Dense.Matrix.Generic.sum x;;

正如我们所看到的，无论身份矩阵中包含什么类型的数字，我们总是将其传递给 Dense.Matrix.Generic.sum 函数。类似地，我们可以对其他模块（Dense.Ndarray、Sparse.Matrix 等）和其他函数（add、mul、neg 等）执行相同的操作。

同时，每个子模块还包含相同的一组函数，例如：

Dense.Matrix.S.(eye 5 |> sum);;

模块快捷方式

实际上，我们经常使用双精度数值，因此 Owl 提供了一些双精度浮点数数据结构的快捷方式：

• Arr 等同于双精度实数 Dense.Ndarray.D；
• Mat 等同于双精度实数 Dense.Matrix.D；

使用这些快捷模块，您不再需要传递类型信息。以下是一些示例。

Arr.zeros [|5|];;        (* 等同于 Dense.Ndarray.D.zeros [|5|] *)
  Mat.zeros 5 5;;          (* 等同于 Dense.Matrix.D.zeros 5 5 *)

除了创建函数之外，还有更多示例。

Mat.load "data.mat";;    (* 等同于 Dense.Matrix.D.load "data.mat" *)
  Mat.of_array 5 5 x;;     (* 等同于 Dense.Matrix.D.of_array 5 5 x *)
  Mat.linspace 0. 9. 10;;  (* 等同于 Dense.Matrix.D.linspace 0. 9. 10 *)

如果您实际上更常使用其他数字类型，比如复数，您当然可以根据需要为相应的 S、D、C 和 Z 模块创建自己的别名。

类型转换

正如我之前提到的，有四种基本数字类型。因此，您可以使用 Generic 模块中的 cast_* 函数将一个值从一种类型转换为另一种类型。

• Generic.cast_s2d: 从 float32 转换为 float64；
• Generic.cast_d2s: 从 float64 转换为 float32；
• Generic.cast_c2z: 从 complex32 转换为 complex64；
• Generic.cast_z2c: 从 complex64 转换为 complex32；
• Generic.cast_s2c: 从 float32 转换为 complex32；
• Generic.cast_d2z: 从 float64 转换为 complex64；
• Generic.cast_s2z: 从 float32 转换为 complex64；
• Generic.cast_d2c: 从 float64 转换为 complex32；

实际上，所有这些函数都依赖于以下 cast 函数。

    val cast : ('a, 'b) kind -> ('c, 'd) t -> ('a, 'b) t
  

第一个参数指定了转换类型。如果源类型和转换类型相同，则 cast 函数简单地复制传递的值。

let x = Arr.uniform [|5;5|]                     (* 在 float64 中创建 *);;
  >val x : Arr.arr =
  >
  >         C0        C1        C2       C3       C4
  >R0 0.648406  0.616945  0.828173 0.579604 0.212017
  >R1 0.960002 0.0563993  0.219521 0.855164 0.414024
  >R2 0.526179  0.532062 0.0640247 0.786426 0.956565
  >R3 0.810557  0.476031  0.516506  0.11439 0.964041
  >R4 0.981665  0.446936  0.276383 0.414747 0.174775
  

现在让我们将 x 从 float64 转换为 complex32。

let y = Dense.Ndarray.Generic.cast Complex32 x  (* 转换为 complex32 *);;
  >val y : (Complex.t, complex32_elt) Dense.Ndarray.Generic.t =
  >
  >               C0              C1              C2             C3             C4
  >R0 (0.648406, 0i)  (0.616945, 0i)  (0.828173, 0i) (0.579604, 0i) (0.212017, 0i)
  >R1 (0.960002, 0i) (0.0563993, 0i)  (0.219521, 0i) (0.855164, 0i) (0.414024, 0i)
  >R2 (0.526179, 0i)  (0.532062, 0i) (0.0640247, 0i) (0.786426, 0i) (0.956565, 0i)
  >R3 (0.810557, 0i)  (0.476031, 0i)  (0.516506, 0i)  (0.11439, 0i) (0.964041, 0i)
  >R4 (0.981664, 0i)  (0.446936, 0i)  (0.276383, 0i) (0.414747, 0i) (0.174775, 0i)
  

要了解每个模块中提供的函数的更多信息，请阅读 Generic 模块的相应接口文件。 Generic 模块包含文档。

• Dense.Ndarray.Generic
• Dense.Matrix.Generic
• Sparse.Ndarray.Generic
• Sparse.Matrix.Generic