User Tools

Site Tools


arraysmatrix


snippet.juliarepl
julia> pkgchk.( [ "julia" => v"1.0.3" ] );

Matrix (2-Dimensional Arrays)

Creating A Matrix from Values

A 2-D array can be created by specifying all elements in blank separated columns and semicolon (not comma!) separated rows:

snippet.juliarepl
julia> x= [ 1 2 3 4 ; 5 6 7 8 ]
2×4 Array{Int64,2}:
 1  2  3  4
 5  6  7  8
  • A newline can be used instead of the semi-colon

WARNING Watch the spaces: x= [ 1 -2 ] is a 1×2 array (matrix). In contrast, x= [ 1-2 ] is a 1-element vector.

A matrix is fundamentally not a vector of vectors. It is a rectangular object. Compare to a vector of vectors:

snippet.juliarepl
julia> x= [ 1, 2, 3 ] ; y= [ 10, 11, 12 ];

julia> [ x, y ]
2-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [10, 11, 12]

Cat (hcat/vcat/hvcat) and Bracket Alternatives

Vectors are created with comma-separated lists. Matrices are created with space and semicolon separated lists.

The parenthesis notation above are “just” shortcuts.

snippet.juliarepl
julia> [ 1 2 3 ] == hcat(1,2,3)
true

julia> [ 1;2;3 ] == vcat(1,2,3)
true

julia> hvcat( (3,3,3), 1,2,3,4,50,6,7,8,9 )
3×3 Array{Int64,2}:
 1  2  3
 4 50  6
 7  8  9

julia> x= collect( 1:5 );  y= collect(11:15); hcat( x, y )	## place multiple vectors into matrix
5×2 Array{Int64,2}:
 1  11
 2  12
 3  13
 4  14
 5  15

Printing A Matrix

snippet.juliarepl
julia> x= [1 2 3 ; 4 5 6];

julia> println(x)
[1 2 3; 4 5 6]

julia> display(x)
2×3 Array{Int64,2}:
 1  2  3
 4  5  6

Initializing a Matrix with Same Values

fill() can initialize arbitrarily-dimensional arrays:

snippet.juliarepl
julia> fill( NaN , (3,2) )	# or use fill( NaN, 3, 2 )
3×2 Array{Float64,2}:
 NaN  NaN
 NaN  NaN
 NaN  NaN

There are type-specialized versions for zeros(Float64, 3, 3) and ones(Bool, 3, 3).

You could also create uninitialized matrices. Even if you immediately assign to your matrix afterwards, such uninitialized variables are rarely a good idea, and the cost of initialization with nada is trivial. Chances are that sooner or later, you may forget. Please always initialize your variables!

Diagonal Matrices

The identity matrix has special support:

snippet.juliarepl
julia> using LinearAlgebra

julia> Matrix( 2.0I, 3, 3 )		## An 'I' is following '2.0' creates a particular diagonal matrix
3×3 Array{Float64,2}:
 2.0  0.0  0.0
 0.0  2.0  0.0
 0.0  0.0  2.0

julia> Matrix( I, 3, 3 )		## I is shortcut for 1.0I
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> Matrix( 10.0I, 3, 3 ) + I	## I is an identity matrix that reshapes itself to whatever is needed
3×3 Array{Float64,2}:
 11.0   0.0   0.0
  0.0  11.0   0.0
  0.0   0.0  11.0

For arbitrary diagonal matrices, you can use

snippet.juliarepl
julia> using LinearAlgebra

julia> diagm( 0 => 1:4 )
4×4 Array{Int64,2}:
 1  0  0  0
 0  2  0  0
 0  0  3  0
 0  0  0  4

Pushing a Vector onto an Array

Remember: A matrix is not like a vector of vectors.

Pushing Vectors onto Vectors

snippet.juliarepl
julia> mtype= Vector{ Vector{Float64} }
Array{Array{Float64,1},1}

julia> M= mtype(undef,0)	## oddly, 'undef' is now required and nothing else is accepted
0-element Array{Array{Float64,1},1}

julia> M= mtype(undef, 2)	## an object of type M with two undef items
2-element Array{Array{Float64,1},1}:
 #undef
 #undef

julia> push!( M, [2.0, 3.0] )	## add a vector to the vector (stack)
3-element Array{Array{Float64,1},1}:
 #undef
 #undef
    [2.0, 3.0]

Side-Append Column Vectors to Matrix

snippet.juliarepl
julia> Matrix{Float64}(undef,2,0)			## Method 1: direct empty 2×0 array: 2 rows, 0 columns
2×0 Array{Float64,2}

julia> M= reshape( Vector{Float64}([]), 2, 0 )       ## Method 2: via reshape of Vector
2×0 Array{Float64,2}

julia> [M [ 1.0, 2.0 ]]                              ## space means hcat. hcat= cat(2,...) , i.e., dimension 2
2×1 Array{Float64,2}:
 1.0
 2.0

julia> [ [M [1.0, 2.0]]  [2.0, 3.0]]                            ## two 2 vectors
2×2 Array{Float64,2}:
 1.0  2.0
 2.0  3.0

julia> M= hcat(M, [ 3.0, 4.0 ] ); M=hcat(M, [ 5.0, 6.0 ] )      ## two 2 vectors
2×2 Array{Float64,2}:
 3.0  5.0
 4.0  6.0

Bottom-Append Row Vectors to Matrix

snippet.juliarepl
julia> M= Matrix{Float64}(undef,0,2)                         ## empty 0×2 array: 2 columns, 0 rows
0×2 Array{Float64,2}

julia> [M;[1.0,2.0]']                                  ## ; means vcat. vcat= cat(1,...) , i.e., dimension 1
1×2 Array{Float64,2}:
 1.0  2.0

julia> vcat( vcat( vcat(M, [1.0,2.0]'), [3.0,4.0]' ), [5.0,6.0]' )  ## push row vector (transposed column vecs)
3×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0
 5.0  6.0

Multidimensional Iterators ("Comprehensions")

A sophisticated “comprehensions” versions of “for” can initialize multidimensional arrays with iterators:

snippet.juliarepl
julia> [r * c for r in 1:5, c in 1:5]
5×5 Array{Int64,2}:
1   2   3   4   5
2   4   6   8  10
3   6   9  12  15
4   8  12  16  20
5  10  15  20  25

Comprehensions are well explained in https://en.wikibooks.org/wiki/Introducing_Julia/Arrays_and_tuples#Creating_arrays_using_range_objects

“Generator Expressions” can omit the surrounding brackets, as in

snippet.juliarepl
julia> sum(1/n^2 for n=1:10)
1.5497677311665408

Matrix Function Over Grid From Two Vectors (for 3D or Contourplots)

Julia makes it easy to create a rectangular grid with values. Here is how the grid will look like:

snippet.juliarepl
julia> x= 1:5; y= 3:9; g(x,y)= (x,y);		## show off the grid that we will generate

julia> g.(x,y')			        	## note element-wise use with transpose
5×7 Array{Tuple{Int64,Int64},2}:
 (1, 3)  (1, 4)  (1, 5)  (1, 6)  (1, 7)  (1, 8)  (1, 9)
 (2, 3)  (2, 4)  (2, 5)  (2, 6)  (2, 7)  (2, 8)  (2, 9)
 (3, 3)  (3, 4)  (3, 5)  (3, 6)  (3, 7)  (3, 8)  (3, 9)
 (4, 3)  (4, 4)  (4, 5)  (4, 6)  (4, 7)  (4, 8)  (4, 9)
 (5, 3)  (5, 4)  (5, 5)  (5, 6)  (5, 7)  (5, 8)  (5, 9)

The prime on the y defines the transpose (or adjoint). Instead of a g() function that returns tuples to illustrate what will be used in each array cell, we can use a simple function:

snippet.juliarepl
julia> x= 1:5; y= 3:9; f(x,y)= x^2*sqrt(y);	## apply a useful function

julia> f.(x,y')
5×7 Array{Float64,2}:
  1.73205   2.0   2.23607   2.44949   2.64575   2.82843   3.0
  6.9282    8.0   8.94427   9.79796  10.583    11.3137   12.0
 15.5885   18.0  20.1246   22.0454   23.8118   25.4558   27.0
 27.7128   32.0  35.7771   39.1918   42.332    45.2548   48.0
 43.3013   50.0  55.9017   61.2372   66.1438   70.7107   75.0

Number of Dimensions, Elements, and Dimensions

snippet.juliarepl
julia> M= [r * c for r in 1:3, c in 1:6]	## a sample two-dimensional comprehension
3×6 Array{Int64,2}:
 1  2  3   4   5   6
 2  4  6   8  10  12
 3  6  9  12  15  18

julia> ndims(M)					## array as two dimensions
2

julia> size(M)				## number of rows (aka nrows, 3) and number of columns (aka ncols, 6)
(3, 6)

julia> size(M,2)			## the second dimension
6

Linear Indexing (Folded Matrix as Vector)

Julia can interpret matrices as long vectors, placed in column-row order. This is called linear indexing:

snippet.juliarepl
julia> M= [r * c for r in 1:3, c in 1:6];	## the same two-dimensional matrix

julia> 



julia> length(M)				## the total number of elements
18

julia> ( M[5], M[8], M[11] )			## 5=(2,2); 8=(2,3); 11=(2,4)
( 4, 6, 8 )

julia> eachindex(M)				## An iterator (OneTo). (could be collect()ed)
Base.OneTo(18)

julia> eachindex(M) == 1:length(M)
true

julia> ( stride(M,1) , stride(M,2) )	## num indexes to advance in vector M to advance in matrix M direction 1 or 2
(1, 3)

Mixing Linear and Non-Linear Access

snippet.juliarepl
julia> M= [ 1 2 3 4 ; 5 6 7 8 ]
2×4 Array{Int64,2}:
 1  2  3  4
 5  6  7  8

julia> M[:]                ## matrix linear access as vector
8-element Array{Int64,1}:
 1
 5
 2
 6
 3
 7
 4
 8

julia> M[:,1]              		## column (same as x[ [1,2] ])
2-element Array{Int64,1}:
 1
 5

julia> M[1,:]              		## row access
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> M[2,3]              		## cell access
7

julia> M[ 2, [3,4] ]       		## arbitrary column and row access
2-element Array{Int64,1}:
 7
 8


julia> M[ :, end–1]        		## `end` is often easier
2-element Array{Int64,1}:
 3
 7

julia> M[ :, [false,true,false,true] ]  ## via Boolean vector
2×2 Array{Int64,2}:
 2  4
 6  8

REPL-Pastable Expression

snippet.juliarepl
julia> M= [r * c for r in 1:3, c in 1:6];        ## the two-dimensional matrix

julia> @show(M);
M = [1 2 3 4 5 6; 2 4 6 8 10 12; 3 6 9 12 15 18]

Converting an X-Dimensional Array into a Different (e.g., Y-Dimensional) Array

To convert an existing 1-D array into a 2-D array, or reshape an existing array, use the reshape function:

snippet.juliarepl
julia> reshape( 1:20 , (4, 5) )              ## sequence into 4×5 matrix
4×5 reshape(::UnitRange{Int64}, 4, 5) with eltype Int64:
 1  5   9  13  17
 2  6  10  14  18
 3  7  11  15  19
 4  8  12  16  20

Usually, the type is ok, but you could also convert( Matrix{Int64}, M ).

WARNING Reshape is an Alias: If you change a value in the reshaped array, the original array will still be affected, too. If this is not desired, work with a copy(v).

snippet.juliarepl
julia> M= [1,2,3,4,5,6];	## the original matrix M

julia> b= reshape(M, 2, 3)	## b is only an alias to M
2×3 Array{Int64,2}:
 1  3  5
 2  4  6

julia> b[1]= 99;		## assigning to b changes M, too!

julia> M			## as you can see here
6-element Array{Int64,1}:
 99
  2
  3
  4
  5
  6

Creating Multidimensional Grids

Use the product function from Iterators:

snippet.juliarepl
julia> using Base.Iterators

julia> arr= collect( product( [1,2,3], [4,5] ) )
3×2 Array{Tuple{Int64,Int64},2}:
 (1, 4)  (1, 5)
 (2, 4)  (2, 5)
 (3, 4)  (3, 5)


julia> reshape( ans, (6,1) )
6×1 Array{Tuple{Int64,Int64},2}:
 (1, 4)
 (2, 4)
 (3, 4)
 (1, 5)
 (2, 5)
 (3, 5)

Tuples were explained in Arrays and Tuples. To untupelize (convert) the resulting object into a matrix, use the reinterpret function (this is a bit tricky due to memory representation formats):

snippet.juliarepl
julia> using Base.Iterators

julia> vt= collect( product( [1,2,3], [4,5] ) )  ## 6-element Array{Tuple{Int64,Int64},1}: vector of 2-tuples
3×2 Array{Tuple{Int64,Int64},2}:
 (1, 4)  (1, 5)
 (2, 4)  (2, 5)
 (3, 4)  (3, 5)

julia> reshape( collect(flatten(vt)), (6,2) )
6×2 Array{Int64,2}:
 1  1
 4  5
 2  2
 4  5
 3  3
 4  5
  • It may not be necessary to materialize the grid with collect, which can save memory.

Learning the Dimension of Higher-Dimensional Arrays

In Julia, rows come first, columns come second. However, this also works for more than two dimensions.

snippet.juliarepl
julia> z= zeros(2,3,4)
2×3×4 Array{Float64,3}:
[:, :, 1] =
 0.0  0.0  0.0
 0.0  0.0  0.0

[:, :, 2] =
 0.0  0.0  0.0
 0.0  0.0  0.0

[:, :, 3] =
 0.0  0.0  0.0
 0.0  0.0  0.0

[:, :, 4] =
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> size( z )
(2, 3, 4)

The second argument to size selects the element from the results vector. Thus, size( zeros(2,3,4), 2 ) == size( zeros(2,3,4) )[2], which here is 3.

snippet.juliarepl
julia> ncol( d ) = size( d, 2 )
ncol (generic function with 1 method)

julia> ncol( zeros(2,3) )
3

Mapping Over Matrices (Applying Operations f(X,Y) -> Z)

For Loops

The standard code is:

snippet.juliarepl
julia> f(x,y) = sqrt( x^2 + y^2 );

julia> xset= 10:11; yset= 12:13; i= 1;

julia> m= fill( NaN, length(xset)*length(yset), 3 )

julia> i
1

julia> for x in xset
	   for y in yset
		global i
		m[i,:]=	[ x, y, f(x,y) ]
		i+=1
	    end#for y#
	end#for x##

julia> m
4×3 Array{Float64,2}:
 10.0  12.0  15.6205
 10.0  13.0  16.4012
 11.0  12.0  16.2788
 11.0  13.0  17.0294

Possible speed enhancements:

  • use Matrix{Float64}(undef, length(xset)*length(yset), 3); to work with unitialized matrix
  • place @inbounds before the for statement tells Julia not to check the bounds. However, this risks crashing Julia and probably works only inside functions.
  • use @simd @inbounds before the for statement to ask for modestly advanced Intel CPU vector features.

Compact Alternatives

snippet.julia
using Base.Iterators

reduce(vcat, [[x, y, f(x, y)]' for (y, x) in Iterators.product(yset, xset)])

or

snippet.julia
hcat([ [z(x,y) for x in xset for y in yset] for z in ((x,y)->x, (x,y)->y, f) ]...)

or

snippet.julia
using StaticArrays

vv = [ @SVector [x, y, f(x,y)]  for x in xset for y in yset ]

reinterpret(Float64, vv, (3, length(vv) ))

or

snippet.julia
using Base.Iterators
reduce(hcat, [[x, y, f(x, y)] for (y, x) in Iterators.product(yset, xset)]) ## 8.531 μs

Broadcasting

dot_postfix_functions operate on scalars inside arrays. To operate on arrays inside arrays (e.g., vectors inside matrices),

snippet.juliarepl
julia> broadcast( / , [1 2 3 4 ; 5 6 7 8], [10,20] )	## shorthand for [1 2 3 4 ; 5 6 7 8 ] ./ [10,20]
2×4 Array{Float64,2}:
 0.1   0.2  0.3   0.4
 0.25  0.3  0.35  0.4

julia> broadcast( / , [1 2 3 4 ; 5 6 7 8], [10,20] )	==  [1 2 3 4 ; 5 6 7 8 ] ./ [10,20]
true

julia> broadcast( / , [10,20], [1 2 3 4 ; 5 6 7 8] )	## shorthand for [10,20] ./ [1 2 3 4 ; 5 6 7 8 ]
2×4 Array{Float64,2}:
 10.0  5.0      3.33333  2.5
  4.0  3.33333  2.85714  2.5

julia> broadcast( / , [10,20], [1 2 3 4 ; 5 6 7 8] )	== [10,20] ./ [1 2 3 4 ; 5 6 7 8 ]
true

The dot function versions are faster than the broadcasts, because the dot functions know how to “fuse” operations. That is, they do not create temporary intermediate matrices—they fuse all the necessary operations into each destination cell.

Removing Sets of Rows Based on Condition (e.g., containing at least one NaN)

snippet.juliarepl
julia> M= collect(reshape( 1.0:18.0, 6, 3 ));  M[2,3]=NaN;  M[4,2]=NaN;  ## rows 4 and 6 are bad

julia> M
6×3 Array{Float64,2}:
 1.0    7.0   13.0
 2.0    8.0  NaN
 3.0    9.0   15.0
 4.0  NaN     16.0
 5.0   11.0   17.0
 6.0   12.0   18.0

julia> goodobs= map( i ->(!any(isnan, M[i,:])), 1:size(M,1) )   ## Method 1
6-element Array{Bool,1}:
  true
 false
  true
 false
  true
  true

julia> M[ goodobs , : ]
4×3 Array{Float64,2}:
 1.0   7.0  13.0
 3.0   9.0  15.0
 5.0  11.0  17.0
 6.0  12.0  18.0


julia> goodobs= filter( i->(all(!isnan, M[i,:])), 1:size(M,1) )		## Method 2
4-element Array{Int64,1}:
 1
 3
 5
 6

julia> M[ goodobs , : ]
4×3 Array{Float64,2}:
 1.0   7.0  13.0
 3.0   9.0  15.0
 5.0  11.0  17.0
 6.0  12.0  18.0

Matrix Overall, Columnwise, and Rowwise Statistics

Some statistics work over matrices. For example,

snippet.juliarepl
julia> using Statistics

julia> M= collect(reshape( 1.0:15.0, 3, 5 ))
3×5 Array{Float64,2}:
 1.0  4.0  7.0  10.0  13.0
 2.0  5.0  8.0  11.0  14.0
 3.0  6.0  9.0  12.0  15.0

julia> ( minimum(M), mean(M), maximum(M) )   ## note: not min() and max(), which operate over their args instead
(1.0, 8.0, 15.0)

julia> mean( M; dims=1 )              ## column-wise row means
1×5 Array{Float64,2}:
 2.0  5.0  8.0  11.0  14.0

julia> std( M; dims=2 )                  ## row-wise column standard deviations
3×1 Array{Float64,2}:
 4.743416490252569
 4.743416490252569
 4.743416490252569

Matrix Math (Linear Algebra)

  • Not shown: inv(A), trace(A), det(A), eigvecs(A), eigvals(A), factorize(A). Many of these, esp factorize, can be clever in exploting aspects of its input matrix.

For matrix multiplication, use the standard *. For element-by-element multiplication, use .*. Powering matrices also works as expected (either matrix with ^, or element-by-element with .^).

snippet.juliarepl
julia> using LinearAlgebra

julia> ones(2,2)^4          ## matrix multiplication, repeated
2×2 Array{Float64,2}:
 8.0  8.0
 8.0  8.0

julia> ones(2,2).^4         ## element multiplicated, repeated
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

Sequential Matrix Multiplications

Consider multiplying four matrices with dimensions (1x10) * (10x10) * (10x10) * (10x1). If you calculate (1x10) * (10x10) and (10x10) * (10x1) first, you need less storage and fewer calculations than if you calculate (10x10) * (10x10) first. Although mathematically identical, the order in which you do calculations matters.

Fortunately, there are two packages that automatically find the speediest way to do sequential matrix multiplications:

Backmatter

Notes

  • Julia has built-in support for sparse arrays, SparseMatrix. Use dropzeros. The implied non-entered value is zero.
  • Julia has built-in support for a variety of special matrices, such as Triangular or Symmetric ones.
  • Julia can easily deal with more than 2-dimensions in arrays (tensors).

References

arraysmatrix.txt · Last modified: 2018/12/28 12:00 (external edit)