Multi-dimensional Arrays¶

Basic Functions¶

Construction and Initialization¶

Many functions for constructing and initializing arrays are provided. In the following list of such functions, calls with a dims... argument can either take a single tuple of dimension sizes or a series of dimension sizes passed as a variable number of arguments.

The syntax [A,B,C,...] constructs a 1-d array (vector) of its arguments.

Concatenation¶

Arrays can be constructed and also concatenated using the following functions:

Scalar values passed to these functions are treated as 1-element arrays.

The concatenation functions are used so often that they have special syntax:

hvcat() concatenates in both dimension 1 (with semicolons) and dimension 2 (with spaces).

Typed array initializers¶

An array with a specific element type can be constructed using the syntax T[A,B,C,...]. This will construct a 1-d array with element type T, initialized to contain elements A, B, C, etc. For example Any[x,y,z] constructs a heterogeneous array that can contain any values.

Concatenation syntax can similarly be prefixed with a type to specify the element type of the result.

julia>[[12][34]]1×4Array{Int64,2}:1234julia>Int8[[12][34]]1×4Array{Int8,2}:1234

Comprehensions¶

Comprehensions provide a general and powerful way to construct arrays. Comprehension syntax is similar to set construction notation in mathematics:

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

The meaning of this form is that F(x,y,...) is evaluated with the variables x, y, etc. taking on each value in their given list of values. Values can be specified as any iterable object, but will commonly be ranges like 1:n or 2:(n-1), or explicit arrays of values like [1.2,3.4,5.7]. The result is an N-d dense array with dimensions that are the concatenation of the dimensions of the variable ranges rx, ry, etc. and each F(x,y,...) evaluation returns a scalar.

The following example computes a weighted average of the current element and its left and right neighbor along a 1-d grid. :

julia>x=rand(8)8-elementArray{Float64,1}:0.8430250.8690520.3651050.6994560.9776530.9949530.410840.809411julia>[0.25*x[i-1]+0.5*x[i]+0.25*x[i+1]fori=2:length(x)-1]6-elementArray{Float64,1}:0.7365590.574680.6854170.9124290.84460.656511

The resulting array type depends on the types of the computed elements. In order to control the type explicitly, a type can be prepended to the comprehension. For example, we could have requested the result in single precision by writing:

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

Generator Expressions¶

Comprehensions can also be written without the enclosing square brackets, producing an object known as a generator. This object can be iterated to produce values on demand, instead of allocating an array and storing them in advance (see Iteration). For example, the following expression sums a series without allocating memory:

julia>sum(1/n^2forn=1:1000)1.6439345666815615

When writing a generator expression with multiple dimensions inside an argument list, parentheses are needed to separate the generator from subsequent arguments:

julia>map(tuple,1/(i+j)fori=1:2,j=1:2,[1:4;])ERROR:syntax:invaliditerationspecification

All comma-separated expressions after for are interpreted as ranges. Adding parentheses lets us add a third argument to map:

julia>map(tuple,(1/(i+j)fori=1:2,j=1:2),[13;24])2×2Array{Tuple{Float64,Int64},2}:(0.5,1)(0.333333,3)(0.333333,2)(0.25,4)

Ranges in generators and comprehensions can depend on previous ranges by writing multiple for keywords:

julia>[(i,j)fori=1:3forj=1:i]6-elementArray{Tuple{Int64,Int64},1}:(1,1)(2,1)(2,2)(3,1)(3,2)(3,3)

In such cases, the result is always 1-d.

Generated values can be filtered using the if keyword:

julia>[(i,j)fori=1:3forj=1:iifi+j==4]2-elementArray{Tuple{Int64,Int64},1}:(2,2)(3,1)

Indexing¶

The general syntax for indexing into an n-dimensional array A is:

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

where each I_k may be:

1. A scalar integer
2. A Range of the form a:b, or a:b:c
3. A : or Colon() to select entire dimensions
4. An arbitrary integer array, including the empty array []
5. A boolean array to select a vector of elements at its true indices

If all the indices are scalars, then the result X is a single element from the array A. Otherwise, X is an array with the same number of dimensions as the sum of the dimensionalities of all the indices.

If all indices are vectors, for example, then the shape of X would be (length(I_1),length(I_2),...,length(I_n)), with location (i_1,i_2,...,i_n) of X containing the value A[I_1[i_1],I_2[i_2],...,I_n[i_n]]. If I_1 is changed to a two-dimensional matrix, then X becomes an n+1-dimensional array of shape (size(I_1,1),size(I_1,2),length(I_2),...,length(I_n)). The matrix adds a dimension. The location (i_1,i_2,i_3,...,i_{n+1}) contains the value at A[I_1[i_1,i_2],I_2[i_3],...,I_n[i_{n+1}]]. All dimensions indexed with scalars are dropped. For example, the result of A[2,I,3] is an array with size size(I). Its ith element is populated by A[2,I[i],3].

Indexing by a boolean array B is effectively the same as indexing by the vector that is returned by find(B). Often referred to as logical indexing, this selects elements at the indices where the values are true, akin to a mask. A logical index must be a vector of the same length as the dimension it indexes into, or it must be the only index provided and match the size and dimensionality of the array it indexes into. It is generally more efficient to use boolean arrays as indices directly instead of first calling find().

Additionally, single elements of a multidimensional array can be indexed as x=A[I], where I is a CartesianIndex. It effectively behaves like an n-tuple of integers spanning multiple dimensions of A. See Iteration below.

As a special part of this syntax, the end keyword may be used to represent the last index of each dimension within the indexing brackets, as determined by the size of the innermost array being indexed. Indexing syntax without the end keyword is equivalent to a call to getindex:

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

Example:

julia>x=reshape(1:16,4,4)4×4Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}}:15913261014371115481216julia>x[2:3,2:end-1]2×2Array{Int64,2}:610711julia>x[map(ispow2,x)]5-elementArray{Int64,1}:124816julia>x[1,[23;41]]2×2Array{Int64,2}:59131

Empty ranges of the form n:n-1 are sometimes used to indicate the inter-index location between n-1 and n. For example, the searchsorted() function uses this convention to indicate the insertion point of a value not found in a sorted array:

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

Assignment¶

The general syntax for assigning values in an n-dimensional array A is:

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

where each I_k may be:

1. A scalar integer
2. A Range of the form a:b, or a:b:c
3. A : or Colon() to select entire dimensions
4. An arbitrary integer array, including the empty array []
5. A boolean array to select elements at its true indices

If X is an array, it must have the same number of elements as the product of the lengths of the indices: prod(length(I_1),length(I_2),...,length(I_n)). The value in location I_1[i_1],I_2[i_2],...,I_n[i_n] of A is overwritten with the value X[i_1,i_2,...,i_n]. If X is not an array, its value is written to all referenced locations of A.

A boolean array used as an index behaves as in getindex(), behaving as though it is first transformed with find().

Index assignment syntax is equivalent to a call to setindex!():

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

Example:

julia>x=collect(reshape(1:9,3,3))3×3Array{Int64,2}:147258369julia>x[1:2,2:3]=-1-1julia>x3×3Array{Int64,2}:1-1-12-1-1369

Iteration¶

The recommended ways to iterate over a whole array are

forainA# Do something with the element aendforiineachindex(A)# Do something with i and/or A[i]end

The first construct is used when you need the value, but not index, of each element. In the second construct, i will be an Int if A is an array type with fast linear indexing; otherwise, it will be a CartesianIndex:

A=rand(4,3)B=view(A,1:3,2:3)julia>foriineachindex(B)@showiendi=Base.IteratorsMD.CartesianIndex_2(1,1)i=Base.IteratorsMD.CartesianIndex_2(2,1)i=Base.IteratorsMD.CartesianIndex_2(3,1)i=Base.IteratorsMD.CartesianIndex_2(1,2)i=Base.IteratorsMD.CartesianIndex_2(2,2)i=Base.IteratorsMD.CartesianIndex_2(3,2)

In contrast with fori=1:length(A), iterating with eachindex provides an efficient way to iterate over any array type.

Array traits¶

If you write a custom AbstractArray type, you can specify that it has fast linear indexing using

Base.linearindexing{T<:MyArray}(::Type{T})=LinearFast()

This setting will cause eachindex iteration over a MyArray to use integers. If you don’t specify this trait, the default value LinearSlow() is used.

Vectorized Operators and Functions¶

The following operators are supported for arrays. The dot version of a binary operator should be used for elementwise operations.

1. Unary arithmetic — -, +, !
2. Binary arithmetic — +, -, *, .*, /, ./, \, .\, ^, .^, div, mod
3. Comparison — .==, .!=, .<, .<=, .>, .>=
4. Unary Boolean or bitwise — ~
5. Binary Boolean or bitwise — &, |, $

Some operators without dots operate elementwise anyway when one argument is a scalar. These operators are *, +, -, and the bitwise operators. The operators / and \ operate elementwise when the denominator is a scalar.

Note that comparisons such as == operate on whole arrays, giving a single boolean answer. Use dot operators for elementwise comparisons.

The following built-in functions are also vectorized, whereby the functions act elementwise:

absabs2anglecbrtairyairyaiairyaiprimeairybiairybiprimeairyprimeacosacoshasinasinhatanatan2atanhacscacschasecasechacotacothcoscospicoshsinsinpisinhtantanhsinccosccsccschsecsechcotcothacosdasindatandasecdacscdacotdcosdsindtandsecdcscdcotdbesselhbesselibesseljbesselj0besselj1besselkbesselybessely0bessely1experferfcerfinverfcinvexp2expm1betadawsondigammaerfcxerfiexponentetazetagammahankelh1hankelh2ceilfloorroundtruncisfiniteisinfisnanlbetalfactlgammaloglog10log1plog2copysignmaxminsignificandsqrthypot

Note that there is a difference between min() and max(), which operate elementwise over multiple array arguments, and minimum() and maximum(), which find the smallest and largest values within an array.

Julia provides the @vectorize_1arg() and @vectorize_2arg() macros to automatically vectorize any function of one or two arguments respectively. Each of these takes two arguments, namely the Type of argument (which is usually chosen to be the most general possible) and the name of the function to vectorize. Here is a simple example:

julia>square(x)=x^2square(genericfunction with1method)julia>@vectorize_1argNumbersquaresquare(genericfunction with2methods)julia>methods(square)# 2 methods for generic function "square":square{T<:Number}(x::AbstractArray{T,N<:Any})atoperators.jl:555square(x)atnone:1julia>square([124;567])2×3Array{Int64,2}:1416253649

Broadcasting¶

It is sometimes useful to perform element-by-element binary operations on arrays of different sizes, such as adding a vector to each column of a matrix. An inefficient way to do this would be to replicate the vector to the size of the matrix:

julia>a=rand(2,1);A=rand(2,3);julia>repmat(a,1,3)+A2×3Array{Float64,2}:1.208131.820681.253871.568511.864011.67846

This is wasteful when dimensions get large, so Julia offers broadcast(), which expands singleton dimensions in array arguments to match the corresponding dimension in the other array without using extra memory, and applies the given function elementwise:

julia>broadcast(+,a,A)2×3Array{Float64,2}:1.208131.820681.253871.568511.864011.67846julia>b=rand(1,2)1×2Array{Float64,2}:0.8675350.00457906julia>broadcast(+,a,b)2×2Array{Float64,2}:1.710560.8476041.736590.873631

Elementwise operators such as .+ and .* perform broadcasting if necessary. There is also a broadcast!() function to specify an explicit destination, and broadcast_getindex() and broadcast_setindex!() that broadcast the indices before indexing. Moreover, f.(args...) is equivalent to broadcast(f,args...), providing a convenient syntax to broadcast any function (Dot Syntax for Vectorizing Functions).

Implementation¶

The base array type in Julia is the abstract type AbstractArray{T,N}. It is parametrized by the number of dimensions N and the element type T. AbstractVector and AbstractMatrix are aliases for the 1-d and 2-d cases. Operations on AbstractArray objects are defined using higher level operators and functions, in a way that is independent of the underlying storage. These operations generally work correctly as a fallback for any specific array implementation.

The AbstractArray type includes anything vaguely array-like, and implementations of it might be quite different from conventional arrays. For example, elements might be computed on request rather than stored. However, any concrete AbstractArray{T,N} type should generally implement at least size(A) (returning an Int tuple), getindex(A,i) and getindex(A,i1,...,iN); mutable arrays should also implement setindex!(). It is recommended that these operations have nearly constant time complexity, or technically Õ(1) complexity, as otherwise some array functions may be unexpectedly slow. Concrete types should also typically provide a similar(A,T=eltype(A),dims=size(A)) method, which is used to allocate a similar array for copy() and other out-of-place operations. No matter how an AbstractArray{T,N} is represented internally, T is the type of object returned by integer indexing (A[1,...,1], when A is not empty) and N should be the length of the tuple returned by size().

DenseArray is an abstract subtype of AbstractArray intended to include all arrays that are laid out at regular offsets in memory, and which can therefore be passed to external C and Fortran functions expecting this memory layout. Subtypes should provide a method stride(A,k) that returns the “stride” of dimension k: increasing the index of dimension k by 1 should increase the index i of getindex(A,i) by stride(A,k). If a pointer conversion method Base.unsafe_convert(Ptr{T},A) is provided, the memory layout should correspond in the same way to these strides.

The Array type is a specific instance of DenseArray where elements are stored in column-major order (see additional notes in Performance Tips). Vector and Matrix are aliases for the 1-d and 2-d cases. Specific operations such as scalar indexing, assignment, and a few other basic storage-specific operations are all that have to be implemented for Array, so that the rest of the array library can be implemented in a generic manner.

SubArray is a specialization of AbstractArray that performs indexing by reference rather than by copying. A SubArray is created with the view() function, which is called the same way as getindex() (with an array and a series of index arguments). The result of view() looks the same as the result of getindex(), except the data is left in place. view() stores the input index vectors in a SubArray object, which can later be used to index the original array indirectly.

StridedVector and StridedMatrix are convenient aliases defined to make it possible for Julia to call a wider range of BLAS and LAPACK functions by passing them either Array or SubArray objects, and thus saving inefficiencies from memory allocation and copying.

The following example computes the QR decomposition of a small section of a larger array, without creating any temporaries, and by calling the appropriate LAPACK function with the right leading dimension size and stride parameters.

julia>a=rand(10,10)10×10Array{Float64,2}:0.5612550.2266780.2033910.308912…0.7503070.2350230.2179640.7189150.5371920.5569460.9962340.6662320.5094230.6607880.4935010.05656220.1183920.4934980.2620480.9406930.2529650.04707790.7369790.2648220.2287870.1614410.8970230.5676410.3439350.323270.7956730.4522420.4688190.6285070.5115280.9355970.9915110.5712970.74485…0.845890.1788340.2844130.1607060.6722520.1331580.655540.3718260.7706280.05312080.3066170.8361260.3011980.02247020.393440.03702050.5360620.8909470.1688770.320020.4861360.0960780.1720480.776720.5077620.5735670.2201240.1658160.2110490.4332770.539476julia>b=view(a,2:2:8,2:2:4)4×2SubArray{Float64,2,Array{Float64,2},Tuple{StepRange{Int64,Int64},StepRange{Int64,Int64}},false}:0.5371920.9962340.7369790.2287870.9915110.744850.8361260.0224702julia>(q,r)=qr(b);julia>q4×2Array{Float64,2}:-0.3388090.78934-0.464815-0.230274-0.6253490.194538-0.527347-0.534856julia>r2×2Array{Float64,2}:-1.58553-0.9215170.00.866567

Compressed Sparse Column (CSC) Storage¶

In Julia, sparse matrices are stored in the Compressed Sparse Column (CSC) format. Julia sparse matrices have the type SparseMatrixCSC{Tv,Ti}, where Tv is the type of the nonzero values, and Ti is the integer type for storing column pointers and row indices.:

type SparseMatrixCSC{Tv,Ti<:Integer}<:AbstractSparseMatrix{Tv,Ti}m::Int# Number of rowsn::Int# Number of columnscolptr::Vector{Ti}# Column i is in colptr[i]:(colptr[i+1]-1)rowval::Vector{Ti}# Row values of nonzerosnzval::Vector{Tv}# Nonzero valuesend

The compressed sparse column storage makes it easy and quick to access the elements in the column of a sparse matrix, whereas accessing the sparse matrix by rows is considerably slower. Operations such as insertion of nonzero values one at a time in the CSC structure tend to be slow. This is because all elements of the sparse matrix that are beyond the point of insertion have to be moved one place over.

All operations on sparse matrices are carefully implemented to exploit the CSC data structure for performance, and to avoid expensive operations.

If you have data in CSC format from a different application or library, and wish to import it in Julia, make sure that you use 1-based indexing. The row indices in every column need to be sorted. If your SparseMatrixCSC object contains unsorted row indices, one quick way to sort them is by doing a double transpose.

In some applications, it is convenient to store explicit zero values in a SparseMatrixCSC. These are accepted by functions in Base (but there is no guarantee that they will be preserved in mutating operations). Such explicitly stored zeros are treated as structural nonzeros by many routines. The nnz() function returns the number of elements explicitly stored in the sparse data structure, including structural nonzeros. In order to count the exact number of actual values that are nonzero, use countnz(), which inspects every stored element of a sparse matrix.

Sparse matrix constructors¶

The simplest way to create sparse matrices is to use functions equivalent to the zeros() and eye() functions that Julia provides for working with dense matrices. To produce sparse matrices instead, you can use the same names with an sp prefix:

julia>spzeros(3,5)3×5sparsematrixwith0Float64nonzeroentriesjulia>speye(3,5)3×5sparsematrixwith3Float64nonzeroentries:[1,1]=1.0[2,2]=1.0[3,3]=1.0

The sparse() function is often a handy way to construct sparse matrices. It takes as its input a vector I of row indices, a vector J of column indices, and a vector V of nonzero values. sparse(I,J,V) constructs a sparse matrix such that 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)5×18sparsematrixwith4Int64nonzeroentries:[1,4]=1[4,7]=2[5,9]=3[3,18]=-5

The inverse of the sparse() function is findn(), which retrieves the inputs used to create the sparse matrix.

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])

Another way to create sparse matrices is to convert a dense matrix into a sparse matrix using the sparse() function:

julia>sparse(eye(5))5×5sparsematrixwith5Float64nonzeroentries:[1,1]=1.0[2,2]=1.0[3,3]=1.0[4,4]=1.0[5,5]=1.0

You can go in the other direction using the full() function. The issparse() function can be used to query if a matrix is sparse.

julia>issparse(speye(5))true

Sparse matrix operations¶

Arithmetic operations on sparse matrices also work as they do on dense matrices. Indexing of, assignment into, and concatenation of sparse matrices work in the same way as dense matrices. Indexing operations, especially assignment, are expensive, when carried out one element at a time. In many cases it may be better to convert the sparse matrix into (I,J,V) format using findnz(), manipulate the non-zeroes or the structure in the dense vectors (I,J,V), and then reconstruct the sparse matrix.

Correspondence of dense and sparse methods¶

The following table gives a correspondence between built-in methods on sparse matrices and their corresponding methods on dense matrix types. In general, methods that generate sparse matrices differ from their dense counterparts in that the resulting matrix follows the same sparsity pattern as a given sparse matrix S, or that the resulting sparse matrix has density d, i.e. each matrix element has a probability d of being non-zero.

Details can be found in the Sparse Vectors and Matrices section of the standard library reference.