How to read this lecture...

Code should execute sequentially if run in a Jupyter notebook

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

a = [10, 20, 30]
3-element Array{Int64,1}:
 10
 20
 30
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

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

size(a)
(3,)
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

eye(3)
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
diagm([2, 4])
2×2 Array{Int64,2}:
 2  0
 0  4
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

Array{Int64, 1} == Vector{Int64}
true
Array{Int64, 2} == Matrix{Int64}
true
Array{Int64, 1} == Matrix{Int64}
false
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()

a = [10, 20, 30, 40]
4-element Array{Int64,1}:
 10
 20
 30
 40
b = reshape(a, 2, 2)
2×2 Array{Int64,2}:
 10  30
 20  40
b
2×2 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:

b[1, 1] = 100  # Continuing the previous example
100
b
2×2 Array{Int64,2}:
 100  30
  20  40
a
4-element Array{Int64,1}:
 100
  20
  30
  40

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

a = [1 2 3 4]  # Two dimensional
1×4 Array{Int64,2}:
 1  2  3  4
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

eye(2)
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0
zeros(3)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

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

x = Array{Float64}(2, 2)
2×2 Array{Float64,2}:
 6.92478e-310  6.92477e-310
 6.92477e-310  0.0

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

ones(2, 2)
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0
fill("foo", 2, 2)
2×2 Array{String,2}:
 "foo"  "foo"
 "foo"  "foo"

Manual Array Definitions

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

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

In two dimensions we can proceed as follows

a = [10 20 30 40]  # Two dimensional, shape is 1 x n
1×4 Array{Int64,2}:
 10  20  30  40
ndims(a)
2
a = [10 20; 30 40]  # 2 x 2
2×2 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

a = [10; 20; 30; 40]
4-element Array{Int64,1}:
 10
 20
 30
 40
ndims(a)
1

Instead transpose the row vector

a = [10 20 30 40]'
4×1 Array{Int64,2}:
 10
 20
 30
 40
ndims(a)
2

Array Indexing

We’ve already seen the basics of array indexing

a = collect(10:10:40)
4-element Array{Int64,1}:
 10
 20
 30
 40
a[end-1]
30
a[1:3]
3-element Array{Int64,1}:
 10
 20
 30

For 2D arrays the index syntax is straightforward

a = randn(2, 2)
2×2 Array{Float64,2}:
 0.665276  -3.12382
 1.49257   -0.756258
a[1, 1]
0.6652764178803451
a[1, :]  # First row
2-element Array{Float64,1}:
  0.665276
 -3.12382
a[:, 1]  # First column
2-element Array{Float64,1}:
 0.665276
 1.49257

Booleans can be used to extract elements

a = randn(2, 2)
2×2 Array{Float64,2}:
 -0.614422   0.0864195
  0.865028  -0.76479
b = [true false; false true]
2×2 Array{Bool,2}:
  true  false
 false   true
a[b]
2-element Array{Float64,1}:
 -0.614422
 -0.76479

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

2-element Array{Float64,1}:
 -0.614422
 -0.76479
a = Array{Float64}(4)
4-element Array{Float64,1}:
 6.92478e-310
 6.92477e-310
 6.92477e-310
 6.92477e-310
a[2:end] = 42
42
a
4-element Array{Float64,1}:
  6.92478e-310
 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

a = ones(3)
3-element Array{Float64,1}:
 1.0
 1.0
 1.0
b = a
3-element Array{Float64,1}:
 1.0
 1.0
 1.0
b[3] = 44
44
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()

a = ones(3)
3-element Array{Float64,1}:
 1.0
 1.0
 1.0
b = copy(a)
3-element Array{Float64,1}:
 1.0
 1.0
 1.0
b[3] = 44
44
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

a = [-1, 0, 1]
3-element Array{Int64,1}:
 -1
  0
  1
length(a)
3
sum(a)
0
mean(a)
0.0
std(a)
1.0
var(a)
1.0
maximum(a)
1
minimum(a)
-1
b = sort(a, rev=true)  # Returns new array, original not modified
3-element Array{Int64,1}:
  1
  0
 -1
b === a  # === tests if arrays are identical (i.e share same memory)
false
b = sort!(a, rev=true)  # Returns *modified original* array
3-element Array{Int64,1}:
  1
  0
 -1
b === a
true

Matrix Algebra

For two dimensional arrays, * means matrix multiplication

a = ones(1, 2)
1×2 Array{Float64,2}:
 1.0  1.0
b = ones(2, 2)
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0
a * b
1×2 Array{Float64,2}:
 2.0  2.0
b * a'
2×1 Array{Float64,2}:
 2.0
 2.0

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

A = [1 2; 2 3]
2×2 Array{Int64,2}:
 1  2
 2  3
B = ones(2, 2)
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0
A \ B
2×2 Array{Float64,2}:
 -1.0  -1.0
  1.0   1.0
inv(A) * B
2×2 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

ones(2) * ones(2)
MethodError: no method matching *(::Array{Float64,1}, ::Array{Float64,1})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...) at operators.jl:138
  *{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},S}(::Union{Base.ReshapedArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2},SubArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:Union{Base.ReshapedArray{T,N,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex,Colon,Int64,Range{Int64}},N}},L}}, ::Union{Base.ReshapedArray{S,1,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray{S,1},SubArray{S,1,A<:Union{Base.ReshapedArray{T,N,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex,Colon,Int64,Range{Int64}},N}},L}}) at linalg/matmul.jl:79
  *(::Base.LinAlg.AbstractTriangular{T,S<:AbstractArray{T,2}}, ::AbstractArray{T,1}) at linalg/triangular.jl:1496
  ...

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

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

b = ones(2, 2)
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0
b * ones(2)
2-element Array{Float64,1}:
 2.0
 2.0
ones(2) * b
DimensionMismatch("A has dimensions (2,1) but B has dimensions (2,2)")

in gemm_wrapper!(::Array{Float64,2}, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,2}) at ./linalg/matmul.jl:309

in *(::Array{Float64,1}, ::Array{Float64,2}) at ./linalg/matmul.jl:86

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

ones(2, 2) * ones(2, 2)   # Matrix multiplication
2×2 Array{Float64,2}:
 2.0  2.0
 2.0  2.0
ones(2, 2) .* ones(2, 2)   # Element by element multiplication
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

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

A = -ones(2, 2)
2×2 Array{Float64,2}:
 -1.0  -1.0
 -1.0  -1.0
A.^2  # Square every element
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

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

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

Scalar multiplication is similar

A = ones(2, 2)
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0
2 * A  # Same as 2 .* A
2×2 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

a = [10, 20, 30]
3-element Array{Int64,1}:
 10
 20
 30
b = [-100, 0, 100]
3-element Array{Int64,1}:
 -100
    0
  100
b .> a
3-element BitArray{1}:
 false
 false
  true
a .== b
3-element BitArray{1}:
 false
 false
 false

We can also do comparisons against scalars with parallel syntax

b
3-element Array{Int64,1}:
 -100
    0
  100
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

a = randn(4)
4-element Array{Float64,1}:
 -1.55242
 -1.0166
 -1.27692
  1.38091
a .< 0
4-element BitArray{1}:
  true
  true
  true
 false
a[a .< 0]
3-element Array{Float64,1}:
 -1.55242
 -1.0166
 -1.27692

Vectorized Functions

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

log(1.0)
0.0

By default, these functions act elementwise on arrays

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

[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

A = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4
det(A)
-2.0
trace(A)
5
eigvals(A)
2-element Array{Float64,1}:
 -0.372281
  5.37228
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}\]

Solutions

Exercise 1

Here’s the iterative approach

function compute_asymptotic_var(A,
                                Sigma,
                                S0=Sigma * Sigma',
                                tolerance=1e-6,
                                maxiter=500)
    V = Sigma * Sigma'
    S = S0
    err = tolerance + 1
    i = 1
    while err > tolerance && i <= maxiter
        next_S = A * S * A' + V
        err = norm(S - next_S)
        S = next_S
        i = i + 1
    end
    return S
end
A =     [0.8 -0.2;
        -0.1 0.7]
Sigma = [0.5 0.4;
         0.4 0.6]

Note that all eigenvalues of \(A\) lie inside the unit disc:

maximum(abs(eigvals(A)))

Let’s compute the asymptotic variance:

0.9
compute_asymptotic_var(A, Sigma)
2×2 Array{Float64,2}:
 0.671228  0.633476
 0.633476  0.858874

Now let’s do the same thing using QuantEcon’s solve_discrete_lyapunov() function and check we get the same result

using QuantEcon
solve_discrete_lyapunov(A, Sigma * Sigma')
2×2 Array{Float64,2}:
 0.671231  0.633474
 0.633474  0.858874