# Vectors, Arrays and Matrices¶

Contents

“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

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

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

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