Linear algebra in Julia: Introduction

2021-01-29

Julia is a relatively new general-purpose programming language (almost 9 years old at the time I write this) that is gaining ground in scientific computation and data science.

This post is the start of an exploration of using Julia to illustrate computations and concepts in a first linear algebra course. The plan is basically to just dive into the linear algebra content and pick up Julia bits (syntax, concepts, etc.) as needed along the way.

Setup

To follow along, you can download and install it on your own machine or get access to it on something like repl.it (account required, the free Starter account should be more than enough for what we'll be doing with Julia; however, note that their version of Julia seems to lag behind, and the free account doesn't give you any privacy). I prefer to run things on my own metal, so it's just a matter of pacman -S julia.

Hello Julia


❯ julia

   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.5.3 (2020-11-09)
 _/ |\__'_|_|_|\__'_|  |
|__/                   |

julia>

Sanity check:


julia> 2+2
4

julia> x = 3
3

julia> x
3

Let's get a 2x2 matrix going:


julia> A = [1 x
            8 -7]
2×2 Array{Int64,2}:
 1   3
 8  -7

Observe that this is called an array and that Julia decided that its entries are integers (fair enough).

The same result can be achieved via


julia> A = [1 x; 8 -7]
2×2 Array{Int64,2}:
 1   3
 8  -7

or


julia> A = [[1, 8] [x, -7]]
2×2 Array{Int64,2}:
 1   3
 8  -7

Note the conspicuous absence of a comma between [1, 8] and [x, -7].

I want a column vector with floating point entries:


julia> v = [1, -6, 9.3, 4]
4-element Array{Float64,1}:
  1.0
 -6.0
  9.3
  4.0

How about complex floats? TODO

Arithmetic

Let's spice things up a bit.


julia> A = [[1, -2, 1] [4, 0, 3]];

julia> B = [[-8, 1, 0] [1, 9, 3]];

julia> A + B
3×2 Array{Int64,2}:
 -7  5
 -1  9
  1  6

julia> C = [[0, 1, -1] [-1, 3, -6] [3, -4, 3]];

julia> 3*C
3×3 Array{Int64,2}:
  0   -3    9
  3    9  -12
 -3  -18    9

julia> 3C
3×3 Array{Int64,2}:
  0   -3    9
  3    9  -12
 -3  -18    9

The latter is kind of nice. However


julia> D = [[8, 0, 1] [-1, 1, 8] [2, -5, 1]];

julia> CD
ERROR: UndefVarError: CD not defined
Stacktrace:
 [1] top-level scope
 [2] run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.3/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288

because Julia takes CD to be the name of a variable, and it has not yet been defined. So we have to do


julia> C*D
3×3 Array{Int64,2}:
  3   23    8
  4  -30  -17
 -5   19   31

Beware the order of multiplication


julia> A*C
ERROR: DimensionMismatch("matrix A has dimensions (3,2), matrix B has dimensions (3,3)")
Stacktrace: (omitted)

julia> C*A
3×2 Array{Int64,2}:
  5   9
 -9  -8
 14   5

Also useful:


julia> transpose(A)
2×3 LinearAlgebra.Transpose{Int64,Array{Int64,2}}:
 1  -2  1
 4   0  3

Up for a stupid trick?


julia> A + 2
ERROR: MethodError: no method matching +(::Array{Int64,2}, ::Int64)
For element-wise addition, use broadcasting with dot syntax: array .+ scalar
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at operators.jl:538
  +(::Complex{Bool}, ::Real) at complex.jl:301
  +(::Missing, ::Number) at missing.jl:115
  ...
Stacktrace: (omitted)

Thank you, Julia, for not following MATLAB's example in miserably failing freshman mathematics. And do I detect a hint of annoyance at MATLAB fandom asking why the operation didn't work and how oh how could one recover this appalling bit of syntactic sugar? Anyway, let's listen to the instruction in the error message:


julia> A .+ 2
3×2 Array{Int64,2}:
 3  6
 0  2
 3  5

Access to entries

Single entries can be accessed for both reading and writing as follows:


julia> A
3×2 Array{Int64,2}:
  1  4
 -2  0
  1  3

julia> A[1, 2]
4

julia> A[1, 2] = 5
5

julia> A
3×2 Array{Int64,2}:
  1  5
 -2  0
  1  3

So indexing is 1-based rather than 0-based.

We can get a whole column or row at a time:


julia> A[:,2]
3-element Array{Int64,1}:
 5
 0
 3

julia> A[1,:]
2-element Array{Int64,1}:
 1
 5

Or maybe you're after a submatrix:


julia> A[2:3,1:2]
2×2 Array{Int64,2}:
 -2  0
  1  3

This is just scratching the surface of Julia's indexing system.

It does bring up the notion of range, so let's play with that:


julia> 1:5
1:5

You don't say! Let's try a little harder:


julia> r = 1:5;

julia> length(r)
5

julia> r[1]
1

julia> r[2]
2

julia> r[3]
3

julia> r[4]
4

julia> r[5]
5

julia> r[6]
ERROR: BoundsError: attempt to access 5-element UnitRange{Int64} at index [6]
Stacktrace: (omitted)

julia> r[end]
5

If you're used to Python, note that in Julia the range contains the endpoint.

Another example:


julia> r = 1:.5:3
1.0:0.5:3.0

julia> length(r)
5

julia> Array(r)
5-element Array{Float64,1}:
 1.0
 1.5
 2.0
 2.5
 3.0

julia> s = 1:.5:3.2
1.0:0.5:3.0

julia> r == s
true

Here the == operator is comparison.

Some matrix creation functions


julia> zeros(3, 2)
3×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0
 0.0  0.0

Bummer, I wanted integers:


julia> zeros(Int, 3, 2)
3×2 Array{Int64,2}:
 0  0
 0  0
 0  0

Better. But let's observe that these are not integers as mathematicians know and love them: they're Int64's, in other words they go from -2^63 to 2^63-1. What comes after that? Let's see:


julia> a = 2^63-1
9223372036854775807

julia> a+1
-9223372036854775808

If you grew up on C or somesuch, this is normal behaviour. If you're coming from Python or the like, you can have actual integers using BigInt


julia> a = BigInt(2^63-1)
9223372036854775807

julia> a+1
9223372036854775808

Let's resume the matrix creation show:


julia> ones(BigInt, 2, 3)
2×3 Array{BigInt,2}:
 1  1  1
 1  1  1

In case you're wondering if there's a function called twelves: no, of course not. But you can do


julia> fill(12, 3, 2)
3×2 Array{Int64,2}:
 12  12
 12  12
 12  12

Identity matrices require pulling in the LinearAlgebra module


julia> Matrix(I, 3, 3)
ERROR: UndefVarError: I not defined
Stacktrace: (omitted)

julia> using LinearAlgebra

julia> Matrix(I, 3, 3)
3×3 Array{Bool,2}:
 1  0  0
 0  1  0
 0  0  1

Seriously, the type is Bool? We can fix this:


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

Another useful construct (also needs LinearAlgebra, but we've already imported that)


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

You can also put smaller matrices together to get larger ones:


julia> [zeros(3, 2) ones(3, 5)]
3×7 Array{Float64,2}:
 0.0  0.0  1.0  1.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0

julia> [zeros(2, 3); ones(5, 3)]
7×3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0