More about matrices

Previously we looked at simple operations with two-dimensional arrays and two-dimensional matrices.

In this section I will also discuss some further topics to cover multi-dimensional arrays, broadcasting and of handling sparse arrays but first I wish to introduce the concept of vectorized and devectorized code.

Vectorized and devectorized code

Consider the following code to add two vectors:

function vecadd1(a,b,c,N)
 for i = 1:N
 c = a + b 
 end
end

function vecadd2(a,b,c,N)
 for i = 1:N, j = 1:length(c)
 c[j] = a[j] + b[j]
 end
end

julia> A = rand(2); B = rand(2); C = zeros(2);
julia> @elapsed vecadd1(A,B,C,100000000)
@elapsed vecadd1(A,B,C,100000000) # => 18.418755286

julia> @elapsed vecadd2(A,B,C,100000000)
@elapsed vecadd2(A,B,C,100000000) # => 0.524002398

Why the difference in timings? The function vecadd1() uses the array plus operation to perform the calculation whereas vecadd2() explicitly loops through the arrays and performs a series of scalar additions. The former is an example of vectorized coding and the latter devectorized, the current situation in Julia is that devectorized code is much quicker than vectorized.

With languages such as R, MATLAB, and Python (using NumPy) vectorized code is faster than devectorized but the reverse is the case in Julia. The reason is that in R (say) vectorization is actually a thin-wrapper around native-C code and since Julia performed is similar to C, calculations which are essentially concerned JUST with array operations will be comparable with those in Julia.

There is little doubt that coding with vector operations is neater and more readable and the designers of Julia are aware of the benefit of improving on timings for vector operations. That it has not been done is tantamount to the difficulty in optimizing code under all circumstances.

Multidimensional arrays

To illustrate the indexing and representation of arrays of higher dimensions, let us generate a three-dimensional array of 64 random numbers laid out in a 4 by 4 by 4 arrangement:

julia> A = rand(4,4,4)
4x4x4 Array{Float64,3}:
[:, :, 1] =
0.522564 0.852847 0.452363 0.444234
0.992522 0.450827 0.885484 0.0693068
0.378972 0.365945 0.757072 0.807745
0.383636 0.383711 0.304271 0.389717

[:, :, 2] =
0.570806 0.912306 0.358262 0.494621
0.810382 0.235757 0.926146 0.915814
0.634989 0.196174 0.773742 0.158593
0.700649 0.843975 0.321075 0.306428

[:, :, 3] =
0.638391 0.606747 0.15706 0.241825
0.492206 0.798426 0.86354 0.715799
0.971428 0.200663 0.00568161 0.0868379
0.936388 0.183021 0.0476718 0.917008

[:, :, 4] =
0.252962 0.432026 0.817504 0.274034
0.164883 0.209135 0.925754 0.876917
0.125772 0.998318 0.593097 0.614772
0.865795 0.204839 0.315774 0.520044

Note

• Use of slice : to display the 3-D matrix
• Can reshape this into a 8x8 2-D matrix
• Values are ordered by the third index, then the second and finally the first

Notice that this can be recast into a two-dimensional array (that is, a matrix) by using the reshape() function:

julia> B = reshape(A,8,8)
8x8 Array{Float64,2}:
0.522564 0.452363 0.570806 0.358262 … 0.15706 0.252962 0.817504
0.992522 0.885484 0.810382 0.926146 0.86354 0.164883 0.925754
0.378972 0.757072 0.634989 0.773742 0.00568161 0.125772 0.593097
0.383636 0.304271 0.700649 0.321075 0.0476718 0.865795 0.315774
0.852847 0.444234 0.912306 0.494621 0.241825 0.432026 0.274034
0.450827 0.0693068 0.235757 0.915814 … 0.715799 0.209135 0.876917
0.365945 0.807745 0.196174 0.158593 0.0868379 0.998318 0.614772
0.383711 0.389717 0.843975 0.306428 0.917008 0.204839 0.520044

Similarly we can reshape the array A into a one-dimensional array:

julia> C = reshape(A,64); typeof(C) # => Array{Float64,1}
julia> transpose(C) ; # we can also write this as C'
1x64 Array{Float64,2}:
0.522564 0.992522 0.378972 0.383636 … 0.876917 0.614772 0.520044

The reshaped array C is of a single dimension but its transpose has two. So the former is a vector and the latter a matrix.

Broadcasting

We saw earlier that it was possible to carry out binary operations on inpidual elements of two arrays by using the operators .+. .* , and so on.

It is sometimes useful to perform operations on arrays of different sizes, such as adding a vector to each column of a matrix. One way to do this is to replicate the vector to the size of the matrix and then apply a .+ operation but this becomes inefficient when the matrix is large.

Julia provides a function 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 element-wise.

The following code generates a 4x3 matrix of Gaussian random numbers, adds 1 to each element and then raises elements of each row to the power 1,2,3,4 respectively:

julia> P = randn(4,3)
4x3 Array{Float64,2}:
 0.859635 -1.05159 1.05167
 0.680754 -1.97133 0.556587
 -0.913965 1.05069 -0.215938
 0.165775 1.72851 -0.884942

julia> Q = [1,2,3,4];
4-element Array{Int64,1}:

julia> broadcast(^,P + 1.0,Q)
4x3 Array{Float64,2}:
 1.85963 -0.0515925 2.05167
 2.82494 0.943481 2.42296
 0.000636823 8.62387 0.482005
 1.84697 55.4247 0.000175252

Sparse matrices

Normal matrices are sometimes referred to as 'dense', which means that there is an entry for cell[i,j]. In cases where most cell values are, say, 0 this is inefficient and it is better to implement a scheme of tuples: (i,j,x) where x is the value referenced by i and j.

These are termed sparse matrices and we can create a sparse matrix by:

S1 = sparse(I, J, X[, m, n, combine])
S2 = sparsevec(I, X[, m, combine])
S3 = sparsevec(D::Dict[, m])

Where S1 of will dimensions m x n and S[I[k], J[k]] = X[k].

If m and n are given they default to max(I) and max(J). The combine function is used to combine duplicates and if not provided, duplicates are added by default.

S2 is a special case where a sparse vector is created and S3 uses an associative array (dictionary) to provide the same thing. The sparse vector is actually an m x 1 size matrix and in the case of S3 row values are keys from the dictionary and the nonzero values are the values from the dictionary (see the section Dictionaries, sets, and others for more information on associative arrays).

Sparse matrices support much of the same set of operations as dense matrices but there are a few special functions which can be applied. For example spzeros(), spones, speye() are the counterparts of zeros(), ones(), and eye() and random number arrays can be generated by sprand() and sprandn():

julia> T = sprand(5,5,0.1)
# The 0.1 means only 10% for the numbers generated will be deemed as nonzero
5x5 sparse matrix with 3 Float64 nonzeros:
 [1, 1] = 0.0942311
 [3, 3] = 0.0165916
 [4, 5] = 0.179855

julia> T*T
5x5 sparse matrix with 2 Float64 non-zeros:
 [1, 1] = 0.00887949
 [3, 3] = 0.000275282

The function full() converts the sparse matrix to a dense one:

julia> S = full(T); typeof(S) # => 5x5 Array{Float64,2}