Vectors, Arrays and Matrices

“Let’s be clear: the work of science has nothing whatever to do with consensus. Consensus is the business of politics. Science, on the contrary, requires only one investigator who happens to be right, which means that he or she has results that are verifiable by reference to the real world. In science consensus is irrelevant. What is relevant is reproducible results.” – Michael Crichton

Overview

In Julia, arrays are the most important data type for working with collections of numerical data

In this lecture we give more details on

  • creating and manipulating Julia arrays
  • fundamental array processing operations
  • basic matrix algebra

Array Basics

Shape and Dimension

We’ve already seen some Julia arrays in action

julia> a = [10, 20, 30]
3-element Array{Int64,1}:
 10
 20
 30

julia> a = ["foo", "bar", 10]
3-element Array{Any,1}:
   "foo"
   "bar"
 10

The REPL tells us that the arrays are of types Array{Int64,1} and Array{Any,1} respectively

Here Int64 and Any are types for the elements inferred by the compiler

We’ll talk more about types later on

The 1 in Array{Int64,1} and Array{Any,1} indicates that the array is one dimensional

This is the default for many Julia functions that create arrays

julia> typeof(randn(100))
Array{Float64,1}

To say that an array is one dimensional is to say that it is flat — neither a row nor a column vector

We can also confirm that a is flat using the size() or ndims() functions

julia> size(a)
(3,)

julia> ndims(a)
1

The syntax (3,) displays a tuple containing one element — the size along the one dimension that exists

Here’s a function that creates a two-dimensional array

julia> eye(3)
3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> diagm([2, 4])
2x2 Array{Int64,2}:
 2  0
 0  4

julia> size(eye(3))
(3,3)

Array vs Vector vs Matrix

In Julia, in addition to arrays you will see the types Vector and Matrix

However, these are just aliases for one- and two-dimensional arrays respectively

julia> Array{Int64, 1} == Vector{Int64}
true

julia> Array{Int64, 2} == Matrix{Int64}
true

julia> Array{Int64, 1} == Matrix{Int64}
false

julia> Array{Int64, 3} == Matrix{Int64}
false

In particular, a Vector in Julia is a flat array

Changing Dimensions

The primary function for changing the dimension of an array is reshape()

julia> a = [10, 20, 30, 40]
4-element Array{Int64,1}:
 10
 20
 30
 40

julia> b = reshape(a, 2, 2)
2x2 Array{Int64,2}:
 10  30
 20  40

julia> b
2x2 Array{Int64,2}:
 10  30
 20  40

Notice that this function returns a “view” on the existing array

This means that changing the data in the new array will modify the data in the old one:

julia> b[1, 1] = 100  # Continuing the previous example
100

julia> b
2x2 Array{Int64,2}:
 100  30
  20  40

julia> a  # First element has changed
4-element Array{Int64,1}:
 100
  20
  30
  40

To collapse an array along one dimension you can use squeeze()

julia> a = [1 2 3 4]  # Two dimensional
1x4 Array{Int64,2}:
 1  2  3  4

julia> squeeze(a, 1)
4-element Array{Int64,1}:
 1
 2
 3
 4

The return value is an Array with the specified dimension “flattened”

Why Flat Arrays?

As we’ve seen, in Julia we have both

  • one-dimensional arrays (i.e., flat arrays)
  • arrays of size (1, n) or (n, 1) that represent row and column vectors respectively

Why do we need both?

On one hand, dimension matters when we come to matrix algebra

  • Multiplying by a row vector is different to multiplication by a column vector

On the other, we use arrays in many settings that don’t involve matrix algebra

In such cases, we don’t care about the distinction between row and column vectors

This is why many Julia functions return flat arrays by default

Creating Arrays

Functions that Return Arrays

We’ve already seen some functions for creating arrays

julia> eye(2)
2x2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0

julia> zeros(3)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

You can create an empty array using the Array() constructor

julia> x = Array(Float64, 2, 2)
2x2 Array{Float64,2}:
 0.0           2.82622e-316
 2.76235e-318  2.82622e-316

The printed values you see here are just garbage values

(the existing contents of the allocated memory slots being interpreted as 64 bit floats)

Other important functions that return arrays are

julia> ones(2, 2)
2x2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0


julia> fill("foo", 2, 2)
2x2 Array{ASCIIString,2}:
 "foo"  "foo"
 "foo"  "foo"

Manual Array Definitions

As we’ve seen, you can create one dimensional arrays from manually specified data like so

julia> a = [10, 20, 30, 40]
4-element Array{Int64,1}:
 10
 20
 30
 40

In two dimensions we can proceed as follows

julia> a = [10 20 30 40]  # Two dimensional, shape is 1 x n
1x4 Array{Int64,2}:
 10  20  30  40

julia> ndims(a)
2

julia> a = [10 20; 30 40]  # 2 x 2
2x2 Array{Int64,2}:
 10  20
 30  40

You might then assume that a = [10; 20; 30; 40] creates a two dimensional column vector but unfortunately this isn’t the case

julia> a = [10; 20; 30; 40]
4-element Array{Int64,1}:
 10
 20
 30
 40

julia> ndims(a)
1

Instead transpose the row vector

julia> a = [10 20 30 40]'
4x1 Array{Int64,2}:
 10
 20
 30
 40

julia> ndims(a)
2

Array Indexing

We’ve already seen the basics of array indexing

julia> a = collect(10:10:40)
4-element Array{Int64,1}:
 10
 20
 30
 40

julia> a[end-1]
30

julia> a[1:3]
3-element Array{Int64,1}:
 10
 20
 30

For 2D arrays the index syntax is straightforward

julia> a = randn(2, 2)
2x2 Array{Float64,2}:
 1.37556  0.924224
 1.52899  0.815694

julia> a[1, 1]
1.375559922478634

julia> a[1, :]  # First row
1x2 Array{Float64,2}:
 1.37556  0.924224

julia> a[:, 1]  # First column
2-element Array{Float64,1}:
 1.37556
 1.52899

Booleans can be used to extract elements

julia> a = randn(2, 2)
2x2 Array{Float64,2}:
 -0.121311  0.654559
 -0.297859  0.89208

julia> b = [true false; false true]
2x2 Array{Bool,2}:
  true  false
 false   true

julia> a[b]
2-element Array{Float64,1}:
 -0.121311
  0.89208

This is useful for conditional extraction, as we’ll see below

An aside: some or all elements of an array can be set equal to one number using slice notation

julia> a = Array(Float64, 4)
4-element Array{Float64,1}:
 1.30822e-282
 1.2732e-313
 4.48229e-316
 1.30824e-282

julia> a[2:end] = 42
42

julia> a
4-element Array{Float64,1}:
  1.30822e-282
 42.0
 42.0
 42.0

Passing Arrays

As in Python, all arrays are passed by reference

What this means is that if a is an array and we set b = a then a and b point to exactly the same data

Hence any change in b is reflected in a

julia> a = ones(3)
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

julia> b = a
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

julia> b[3] = 44
44

julia> a
3-element Array{Float64,1}:
  1.0
  1.0
 44.0

If you are a MATLAB programmer perhaps you are recoiling in horror at this idea

But this is actually the more sensible default – after all, it’s very inefficient to copy arrays unnecessarily

If you do need an actual copy in Julia, just use copy()

julia> a = ones(3)
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

julia> b = copy(a)
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

julia> b[3] = 44
44

julia> a
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

Operations on Arrays

Array Methods

Julia provides standard functions for acting on arrays, some of which we’ve already seen

julia> a = [-1, 0, 1]
3-element Array{Int64,1}:
 -1
  0
  1

julia> length(a)
3

julia> sum(a)
0

julia> mean(a)
0.0

julia> std(a)
1.0

julia> var(a)
1.0

julia> maximum(a)
1

julia> minimum(a)
-1

julia> b = sort(a, rev=true)  # Returns new array, original not modified
3-element Array{Int64,1}:
  1
  0
 -1

julia> b === a  # === tests if arrays are identical (i.e share same memory)
false

julia> b = sort!(a, rev=true)  # Returns *modified original* array
3-element Array{Int64,1}:
  1
  0
 -1

julia> b === a
true

Matrix Algebra

For two dimensional arrays, * means matrix multiplication

julia> a = ones(1, 2)
1x2 Array{Float64,2}:
 1.0  1.0

julia> b = ones(2, 2)
2x2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

julia> a * b
1x2 Array{Float64,2}:
 2.0  2.0

julia> b * a'
2x1 Array{Float64,2}:
 2.0
 2.0

To solve the linear system A X = B for X use A \ B

julia> A = [1 2; 2 3]
2x2 Array{Int64,2}:
 1  2
 2  3

julia> B = ones(2, 2)
2x2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

julia> A \ B
2x2 Array{Float64,2}:
 -1.0  -1.0
  1.0   1.0

julia> inv(A) * B
2x2 Array{Float64,2}:
 -1.0  -1.0
  1.0   1.0

Although the last two operations give the same result, the first one is numerically more stable and should be preferred in most cases

Multiplying two one dimensional vectors gives an error — which is reasonable since the meaning is ambiguous

julia> ones(2) * ones(2)
ERROR: `*` has no method matching *(::Array{Float64,1}, ::Array{Float64,1})

If you want an inner product in this setting use dot()

julia> dot(ones(2), ones(2))
2.0

Matrix multiplication using one dimensional vectors is a bit inconsistent — pre-multiplication by the matrix is OK, but post-multiplication gives an error

julia> b = ones(2, 2)
2x2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

julia> b * ones(2)
2-element Array{Float64,1}:
 2.0
 2.0

julia> ones(2) * b
ERROR: DimensionMismatch("*")
 in gemm_wrapper! at linalg/matmul.jl:275
 in * at linalg/matmul.jl:74

It’s probably best to give your vectors dimension before you multiply them against matrices

Elementwise Operations

Algebraic Operations

Suppose that we wish to multiply every element of matrix A with the corresponding element of matrix B

In that case we need to replace * (matrix multiplication) with .* (elementwise multiplication)

For example, compare

julia> ones(2, 2) * ones(2, 2)   # Matrix multiplication
2x2 Array{Float64,2}:
 2.0  2.0
 2.0  2.0

julia> ones(2, 2) .* ones(2, 2)   # Element by element multiplication
2x2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

This is a general principle: .x means apply operator x elementwise

julia> A = -ones(2, 2)
2x2 Array{Float64,2}:
 -1.0  -1.0
 -1.0  -1.0

julia> A.^2  # Square every element
2x2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

However in practice some operations are unambiguous and hence the . can be omitted

julia> ones(2, 2) + ones(2, 2)  # Same as ones(2, 2) .+ ones(2, 2)
2x2 Array{Float64,2}:
 2.0  2.0
 2.0  2.0

Scalar multiplication is similar

julia> A = ones(2, 2)
2x2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

julia> 2 * A  # Same as 2 .* A
2x2 Array{Float64,2}:
 2.0  2.0
 2.0  2.0

In fact you can omit the * altogether and just write 2A

Elementwise Comparisons

Elementwise comparisons also use the .x style notation

julia> a = [10, 20, 30]
3-element Array{Int64,1}:
 10
 20
 30

julia> b = [-100, 0, 100]
3-element Array{Int64,1}:
 -100
    0
  100

julia> b .> a
3-element BitArray{1}:
 false
 false
  true

julia> a .== b
3-element BitArray{1}:
 false
 false
 false

We can also do comparisons against scalars with parallel syntax

julia> b
3-element Array{Int64,1}:
 -100
    0
  100

julia> b .> 1
3-element BitArray{1}:
 false
 false
  true

This is particularly useful for conditional extraction — extracting the elements of an array that satisfy a condition

julia> a = randn(4)
4-element Array{Float64,1}:
  0.0636526
  0.933701
 -0.734085
  0.531825

julia> a .< 0
4-element BitArray{1}:
 false
 false
  true
 false

julia> a[a .< 0]
1-element Array{Float64,1}:
 -0.734085

Vectorized Functions

Julia provides standard mathematical functions such as log, exp, sin, etc.

julia> log(1.0)
0.0

By default, these functions act elementwise on arrays

julia> log(ones(4))
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0

Functions that act elementwise on arrays in this manner are called vectorized functions

Note that we can get the same result as with a comprehension or more explicit loop

julia> [log(x) for x in ones(4)]
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0

In Julia loops are typically fast and hence the need for vectorized functions is less intense than for some other high level languages

Nonetheless the syntax is convenient

Linear Algebra

Julia provides some a great deal of additional functionality related to linear operations

julia> A = [1 2; 3 4]
2x2 Array{Int64,2}:
 1  2
 3  4

julia> det(A)
-2.0

julia> trace(A)
5

julia> eigvals(A)
2-element Array{Float64,1}:
 -0.372281
  5.37228

julia> rank(A)
2

For more details see the linear algebra section of the standard library

Exercises

Exercise 1

This exercise is on some matrix operations that arise in certain problems, including when dealing with linear stochastic difference equations

If you aren’t familiar with all the terminology don’t be concerned — you can skim read the background discussion and focus purely on the matrix exercise

With that said, consider the stochastic difference equation

(1)\[X_{t+1} = A X_t + b + \Sigma W_{t+1}\]

Here

  • \(X_t, b\) and \(X_{t+1}\) ar \(n \times 1\)
  • \(A\) is \(n \times n\)
  • \(\Sigma\) is \(n \times k\)
  • \(W_t\) is \(k \times 1\) and \(\{W_t\}\) is iid with zero mean and variance-covariance matrix equal to the identity matrix

Let \(S_t\) denote the \(n \times n\) variance-covariance matrix of \(X_t\)

Using the rules for computing variances in matrix expressions, it can be shown from (1) that \(\{S_t\}\) obeys

(2)\[S_{t+1} = A S_t A' + \Sigma \Sigma'\]

It can be shown that, provided all eigenvalues of \(A\) lie within the unit circle, the sequence \(\{S_t\}\) converges to a unique limit \(S\)

This is the unconditional variance or asymptotic variance of the stochastic difference equation

As an exercise, try writing a simple function that solves for the limit \(S\) by iterating on (2) given \(A\) and \(\Sigma\)

To test your solution, observe that the limit \(S\) is a solution to the matrix equation

(3)\[S = A S A' + Q \quad \text{where} \quad Q := \Sigma \Sigma'\]

This kind of equation is known as a discrete time Lyapunov equation

The QuantEcon package provides a function called solve_discrete_lyapunov that implements a fast “doubling” algorithm to solve this equation

Test your iterative method against solve_discrete_lyapunov using matrices

\[\begin{split}A = \begin{bmatrix} 0.8 & -0.2 \\ -0.1 & 0.7 \end{bmatrix} \qquad \Sigma = \begin{bmatrix} 0.5 & 0.4 \\ 0.4 & 0.6 \end{bmatrix}\end{split}\]