Julia 入门教程 Julia 多维数组

多维数组

数组是一个存在多维网格中的对象集合。通常，数组包含的对象的类型为 Any。对大多数计算而言，数组对象一般更具体为 Float64或 Int32。

因为性能的原因，Julia 不希望把程序写成向量化的形式。

在 Julia 中，通过引用将参数传递给函数。Julia 的库函数不会修改传递给它的输入。用户写代码时，如果要想做类似的功能，要注意先把输入复制一份儿。

数组

基础函数

构造和初始化

下列函数中调用的 dims...参数，既可以是维度的单多元组，也可以是维度作为可变参数时的一组值。

连接

使用下列函数，可在任意维度连接数组：

传递给这些函数的标量值被视为一元阵列。

级联功能非常常用，所以为它们设计了特殊的语法：

hvcat可以实现一维上的（使用分号间隔）或二维上的（使用空格间隔）的级联。

Comprehensions

Comprehensions 用于构造数组。它的语法类似于数学中的集合标记法：

A = [ F(x,y,...) for x=rx, y=ry, ... ]

F(x,y,...)根据变量 x, y等来求值。这些变量的值可以是任何迭代对象，但大多数情况下，都使用类似于 1:n或 2:(n-1)的范围对象，或显式指明为类似 [1.2, 3.4, 5.7]的数组。它的结果是 N 维稠密数组。

下例计算在维度 1 上，当前元素及左右邻居元素的加权平均数：

julia> const x = rand(8)
    8-element Array{Float64,1}:
     0.843025
     0.869052
     0.365105
     0.699456
     0.977653
     0.994953
     0.41084
     0.809411
    julia> [ 0.25*x[i-1] + 0.5*x[i] + 0.25*x[i+1] for i=2:length(x)-1 ]
    6-element Array{Float64,1}:
     0.736559
     0.57468
     0.685417
     0.912429
     0.8446
     0.656511

注解： 上例中，x被声明为常量，因为对于非常量的全局变量，Julia 的类型推断不怎么样。

可在 comprehension 之前显式指明它的类型。如要避免在前例中声明 x为常量，但仍要确保结果类型为 Float64，应这样写：

Float64[ 0.25*x[i-1] + 0.5*x[i] + 0.25*x[i+1] for i=2:length(x)-1 ]

使用花括号来替代方括号，可以将它简写为 Any类型的数组：

julia> { i/2 for i = 1:3 }
    3-element Array{Any,1}:
     0.5
     1.0
     1.5

索引

索引 n 维数组 A 的通用语法为：

X = A[I_1, I_2, ..., I_n]

其中 I_k 可以是：

1. 标量
2. 满足 :, a:b, 或 a:b:c格式的 Range对象
3. 任意整数向量，包括空向量 []
4. 布尔值向量

结果 X 的维度通常为 (length(I_1), length(I_2), ..., length(I_n))，且 X 的索引 (i_1, i_2, ..., i_n)处的值为 A[I_1[i_1], I_2[i_2], ..., I_n[i_n]]。缀在后面的标量索引的维度信息被舍弃。如，A[I, 1]的维度为 (length(I),)。布尔值向量先由 find函数进行转换。由布尔值向量索引的维度长度，是向量中 true值的个数。

索引语法与调用 getindex等价：

X = getindex(A, I_1, I_2, ..., I_n)

例如：

julia> x = reshape(1:16, 4, 4)
    4x4 Array{Int64,2}:
     1  5   9  13
     2  6  10  14
     3  7  11  15
     4  8  12  16
    julia> x[2:3, 2:end-1]
    2x2 Array{Int64,2}:
     6  10
     7  11

n:n-1形式的空范围有时用来表示相互索引位置在 n-1和 n之间。例如，在 searchsorted函数使用本习惯指出插入点的值不在排序后的数组中：

julia> a = [1,2,5,6,7];
    julia> searchsorted(a, 3)
    3:2

赋值

给 n 维数组 A 赋值的通用语法为：

A[I_1, I_2, ..., I_n] = X

其中 I_k 可能是：

1. 标量
2. 满足 :, a:b, 或 a:b:c格式的 Range对象
3. 任意整数向量，包括空向量 []
4. 布尔值向量

如果 X是一个数组，它的维度应为 (length(I_1), length(I_2), ..., length(I_n))，且 A在 i_1, i_2, ..., i_n处的值被覆写为 X[I_1[i_1], I_2[i_2], ..., I_n[i_n]]。如果 X不是数组，它的值被写进所有 A被引用的地方。

用于索引的布尔值向量与 getindex中一样（先由 find函数进行转换）。

索引赋值语法等价于调用 setindex!：

setindex!(A, X, I_1, I_2, ..., I_n)

例如：

julia> x = reshape(1:9, 3, 3)
    3x3 Array{Int64,2}:
     1  4  7
     2  5  8
     3  6  9
    julia> x[1:2, 2:3] = -1
    -1
    julia> x
    3x3 Array{Int64,2}:
     1  -1  -1
     2  -1  -1
     3   6   9

向量化的运算符和函数

数组支持下列运算符。逐元素进行的运算，应使用带“点”（逐元素）版本的二元运算符。

1. 一元： -, +, !
2. 二元： +, -, *, .*, /, ./, \, .\, ^, .^, div, mod
3. 比较： .==, .!=, .<, .<=, .>, .>=
4. 一元布尔值或位运算： ~
5. 二元布尔值或位运算： &, |, $

一些没有“点”（逐元素）操作运算符当一个参数是一个标量时会被使用。这些运算符有 *, /, \和按位运算符。

请注意，像 ==操作这样的比较运算是操作在整个阵列上的，它会给出一个布尔返回值。逐位的比较使用点操作符。

下列内置的函数也都是向量化的, 即函数是逐元素版本的：

abs abs2 angle cbrt
    airy airyai airyaiprime airybi airybiprime airyprime
    acos acosh asin asinh atan atan2 atanh
    acsc acsch asec asech acot acoth
    cos  cospi cosh  sin  sinpi sinh  tan  tanh  sinc  cosc
    csc  csch  sec  sech  cot  coth
    acosd asind atand asecd acscd acotd
    cosd  sind  tand  secd  cscd  cotd
    besselh besseli besselj besselj0 besselj1 besselk bessely bessely0 bessely1
    exp  erf  erfc  erfinv erfcinv exp2  expm1
    beta dawson digamma erfcx erfi
    exponent eta zeta gamma
    hankelh1 hankelh2
     ceil  floor  round  trunc
    iceil ifloor iround itrunc
    isfinite isinf isnan
    lbeta lfact lgamma
    log log10 log1p log2
    copysign max min significand
    sqrt hypot

注意 minmax和 minimummaximum之间的区别，前者是对多个数组操作，找出各数组对应的的元素中的最大最小，后者是作用在一个数组上找出该数组的最大最小值。

Julia 提供了 @vectorize_1arg和 @vectorize_2arg两个宏，分别用来向量化任意的单参数或两个参数的函数。每个宏都接收两个参数，即函数参数的类型和函数名。例如：

julia> square(x) = x^2
    square (generic function with 1 method)
    julia> @vectorize_1arg Number square
    square (generic function with 4 methods)
    julia> methods(square)
    # 4 methods for generic function "square":
    square{T<:Number}(::AbstractArray{T<:Number,1}) at operators.jl:359
    square{T<:Number}(::AbstractArray{T<:Number,2}) at operators.jl:360
    square{T<:Number}(::AbstractArray{T<:Number,N}) at operators.jl:362
    square(x) at none:1
    julia> square([1 2 4; 5 6 7])
    2x3 Array{Int64,2}:
      1   4  16
     25  36  49

Broadcasting

有时要对不同维度的数组进行逐元素的二元运算，如将向量加到矩阵的每一列。低效的方法是，把向量复制成同维度的矩阵：

julia> a = rand(2,1); A = rand(2,3);
    julia> repmat(a,1,3)+A
    2x3 Array{Float64,2}:
     1.20813  1.82068  1.25387
     1.56851  1.86401  1.67846

维度很大时，效率会很低。Julia 提供 broadcast函数，它将数组参数的维度进行扩展，使其匹配另一个数组的对应维度，且不需要额外内存，最后再逐元素调用指定的二元函数：

julia> broadcast(+, a, A)
    2x3 Array{Float64,2}:
     1.20813  1.82068  1.25387
     1.56851  1.86401  1.67846
    julia> b = rand(1,2)
    1x2 Array{Float64,2}:
     0.867535  0.00457906
    julia> broadcast(+, a, b)
    2x2 Array{Float64,2}:
     1.71056  0.847604
     1.73659  0.873631

逐元素的运算符，如 .+和 .*将会在必要时进行 broadcasting 。还提供了 broadcast!函数，可以明确指明目的，而 broadcast_getindex和 broadcast_setindex!函数可以在索引前对索引值做 broadcast。

实现

Julia 的基础数组类型是抽象类型 AbstractArray{T,N}，其中维度为 N，元素类型为 T。AbstractVector和 AbstractMatrix分别是它 1 维和 2 维的别名。

AbstractArray类型包含任何形似数组的类型， 而且它的实现和通常的数组会很不一样。例如，任何具体的 AbstractArray{T，N}至少要有 size(A)(返回 Int多元组)，getindex(A,i)和 getindex(A,i1,...,iN)(返回 T类型的一个元素), 可变的数组要能 setindex！。 这些操作都要求在近乎常数的时间复杂度或 O(1) 复杂度，否则某些数组函数就会特别慢。具体的类型也要提供类似于 similar(A,T=eltype(A),dims=size(A))的方法用来分配一个拷贝。

DenseArray是一个抽象的 AbstractArray类型的亚型，它应该包括在内存的常规偏移上的所有数组，因此可以被传递到外部在此内存布局上的 C 和 Fortran 函数。
亚型应该提供一个方法 stride(A,k)，使之返回“跨越”的维度 k：向给出的维度 k加 1应该使 getindex(A,i)中的 i增加 stride(A,k)。如果提供了一个指针转换方法 convert(Ptr{T},A)，那么内存布局应该以相同的方式对应于这些扩展。

Array{T,N}类型是 DenseArray的特殊实例，它的元素以列序为主序存储（详见代码性能优化）。Vector和 Matrix是分别是它 1 维和 2 维的别名。

SubArray是 AbstractArray的特殊实例，它通过引用而不是复制来进行索引。使用 sub函数来构造 SubArray，它的调用方式与 getindex相同（使用数组和一组索引参数）。sub的结果与 getindex的结果类似，但它的数据仍留在原地。sub在 SubArray对象中保存输入的索引向量，这个向量将被用来间接索引原数组。

StridedVector和 StridedMatrix是为了方便而定义的别名。通过给他们传递 Array或 SubArray对象，可以使 Julia 大范围调用 BLAS 和 LAPACK 函数，提高内存申请和复制的效率。

下面的例子计算大数组中的一个小块的 QR 分解，无需构造临时变量，直接调用合适的 LAPACK 函数。

julia> a = rand(10,10)
    10x10 Array{Float64,2}:
     0.561255   0.226678   0.203391  0.308912   …  0.750307  0.235023   0.217964
     0.718915   0.537192   0.556946  0.996234      0.666232  0.509423   0.660788
     0.493501   0.0565622  0.118392  0.493498      0.262048  0.940693   0.252965
     0.0470779  0.736979   0.264822  0.228787      0.161441  0.897023   0.567641
     0.343935   0.32327    0.795673  0.452242      0.468819  0.628507   0.511528
     0.935597   0.991511   0.571297  0.74485    …  0.84589   0.178834   0.284413
     0.160706   0.672252   0.133158  0.65554       0.371826  0.770628   0.0531208
     0.306617   0.836126   0.301198  0.0224702     0.39344   0.0370205  0.536062
     0.890947   0.168877   0.32002   0.486136      0.096078  0.172048   0.77672
     0.507762   0.573567   0.220124  0.165816      0.211049  0.433277   0.539476
    julia> b = sub(a, 2:2:8,2:2:4)
    4x2 SubArray{Float64,2,Array{Float64,2},(StepRange{Int64,Int64},StepRange{Int64,Int64})}:
     0.537192  0.996234
     0.736979  0.228787
     0.991511  0.74485
     0.836126  0.0224702
    julia> (q,r) = qr(b);
    julia> q
    4x2 Array{Float64,2}:
     -0.338809   0.78934
     -0.464815  -0.230274
     -0.625349   0.194538
     -0.527347  -0.534856
    julia> r
    2x2 Array{Float64,2}:
     -1.58553  -0.921517
      0.0       0.866567

稀疏矩阵

稀疏矩阵是其元素大部分为 0 的矩阵。

列压缩（CSC）存储

Julia 中，稀疏矩阵使用 列压缩（CSC）格式。Julia 稀疏矩阵的类型为 SparseMatrixCSC{Tv,Ti}，其中 Tv是非零元素的类型， Ti是整数类型，存储列指针和行索引：

type SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
        m::Int                  # Number of rows
        n::Int                  # Number of columns
        colptr::Vector{Ti}      # Column i is in colptr[i]:(colptr[i+1]-1)
        rowval::Vector{Ti}      # Row values of nonzeros
        nzval::Vector{Tv}       # Nonzero values
    end

列压缩存储便于按列简单快速地存取稀疏矩阵的元素，但按行存取则较慢。把非零值插入 CSC 结构等运算，都比较慢，这是因为稀疏矩阵中，在所插入元素后面的元素，都要逐一移位。

如果你从其他地方获得的数据是 CSC 格式储存的，想用 Julia 来读取，应确保它的序号从 1 开始索引。每一列中的行索引值应该是排好序的。如果你的 SparseMatrixCSC对象包含未排序的行索引值，对它们进行排序的最快的方法是转置两次。

有时，在 SparseMatrixCSC中存储一些零值，后面的运算比较方便。 Base中允许这种行为（但是不保证在操作中会一直保留这些零值）。这些被存储的零被许多函数认为是非零值。nnz函数返回稀疏数据结构中存储的元素数目，包括被存储的零。要想得到准确的非零元素的数目，请使用 countnz函数，它挨个检查每个元素的值（因此它的时间复杂度不再是常数，而是与元素数目成正比）。

构造稀疏矩阵

稠密矩阵有 zeros和 eye函数，稀疏矩阵对应的函数，在函数名前加 sp前缀即可：

julia> spzeros(3,5)
    3x5 sparse matrix with 0 Float64 entries:
    julia> speye(3,5)
    3x5 sparse matrix with 3 Float64 entries:
            [1, 1]  =  1.0
            [2, 2]  =  1.0
            [3, 3]  =  1.0

sparse函数是比较常用的构造稀疏矩阵的方法。它输入行索引 I，列索引向量 J，以及非零值向量 V。 sparse(I,J,V)构造一个满足 S[I[k], J[k]] = V[k]的稀疏矩阵：

julia> I = [1, 4, 3, 5]; J = [4, 7, 18, 9]; V = [1, 2, -5, 3];
    julia> S = sparse(I,J,V)
    5x18 sparse matrix with 4 Int64 entries:
            [1 ,  4]  =  1
            [4 ,  7]  =  2
            [5 ,  9]  =  3
            [3 , 18]  =  -5

与 sparse相反的函数为 findn，它返回构造稀疏矩阵时的输入：

julia> findn(S)
    ([1,4,5,3],[4,7,9,18])
    julia> findnz(S)
    ([1,4,5,3],[4,7,9,18],[1,2,3,-5])

另一个构造稀疏矩阵的方法是，使用 sparse函数将稠密矩阵转换为稀疏矩阵：

julia> sparse(eye(5))
    5x5 sparse matrix with 5 Float64 entries:
            [1, 1]  =  1.0
            [2, 2]  =  1.0
            [3, 3]  =  1.0
            [4, 4]  =  1.0
            [5, 5]  =  1.0

可以使用 dense或 full函数做逆操作。 issparse函数可用来检查矩阵是否稀疏：

julia> issparse(speye(5))
    true

稀疏矩阵运算

稠密矩阵的算术运算也可以用在稀疏矩阵上。对稀疏矩阵进行赋值运算，是比较费资源的。大多数情况下，建议使用 findnz函数把稀疏矩阵转换为 (I,J,V)格式，在非零数或者稠密向量 (I,J,V)的结构上做运算，最后再重构回稀疏矩阵。

稠密矩阵和稀疏矩阵函数对应关系

接下来的表格列出了内置的稀疏矩阵的函数, 及其对应的稠密矩阵的函数。通常，稀疏矩阵的函数，要么返回与输入稀疏矩阵 S同样的稀疏度，要么返回 d稠密度，例如矩阵的每个元素是非零的概率为 d。

详见可以标准库文档的 stdlib-sparse章节。