Linear algebra in Julia: Introduction
A. Ghitza
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