User Tools

Site Tools


arraysmatrix


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

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

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 ]; z= [ x, y ]
2-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [10, 11, 12]

hcat/vcat/hvcat

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

The parenthesis notation above is 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 be used to initialize arbitrarily-dimensional arrays:

snippet.juliarepl
julia> 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), but not fill().

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

Arbitrary Diagonal

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

The Identity Matrix

snippet.juliarepl
julia> using LinearAlgebra

julia> Matrix( 1.0I, 3, 3 )	## look closely.  An 'I' is following '1.0'
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

Pushing a Vector onto an Array

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]

Appending Vectors onto Matrices

Side-Append Column Vectors

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

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 From Two Vectors (for 3D or Contourplots)

snippet.juliarepl
julia> x= [1:5;]; y= [3:9;]; f(x,y)= x^2*sqrt(y);

julia> f.(x,y')        ## note element-wise use with transpose
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]
3×6 Array{Int64,2}:
 1  2  3   4   5   6
 2  4  6   8  10  12
 3  6  9  12  15  18

julia> length(M)
18

julia> ndims(M)
2

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

julia> size(M,2)
6

julia> eachindex(M)
Base.OneTo(18)

julia> collect(eachindex(M)) == collect([1:length(M);])
true

julia> ( stride(M,1) , stride(M,2) )
(1, 3)

Accessing Rows, Columns, or Cells of a Matrix

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 accessed as vector
8-element Array{Int64,1}:
 1
 5
 2
 6
 3
 7
 4
 8

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

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

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

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


julia> M[ :, end–1]        ## `end` is useful
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

To obtain an expression that you can paste into the REPL, use @show M.

Direct Access to Matrix as Folded Vector

snippet.juliarepl
julia> M= [ 1 2 3 4 ; 5 6 7 8 ]; M[2]       ## internally stored in column-row order
5

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];  b= reshape(M, 2, 3)
2×3 Array{Int64,2}:
 1  3  5
 2  4  6

julia> b[1]= 99; M
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

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 “naive” and good-enough code is:

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

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

julia> m= Matrix{Float64}(undef, 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 x
	end#for##

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

For speed, it can be improved to

snippet.julia
    @inbounds for x in Float64.(xset) 
        @simd for y in Float64.(yset)
            m[:,i] .= ( x, y, f(x,y) )
            i += 1
        end
    end

Note that Float64. to help SIMD. The following is also surprisingly fast:

snippet.julia
	m = Matrix{Float64}(undef, 3, length(xset)*length(yset))
	i = 1
	@inbounds for x in xset, y in yset 
		m[1,i] = x
		m[2,i] = y
        	m[3,i] = f(x,y)
		i += 1
	end

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

[[functions|dot_postfix_functions]] operate on scalars inside arrays.  To operate on arrays inside arrays (e.g., vectors inside matrices),

julia> broadcast(/, [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(/, [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

Removing Sets of Rows Based on Condition (e.g., with 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/11/22 20:47 (external edit)