Alex Ghitza

## 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
``````

### 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:
 top-level scope
 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)
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
``````

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

julia> r
2

julia> r
3

julia> r
4

julia> r
5

julia> r
ERROR: BoundsError: attempt to access 5-element UnitRange{Int64} at index 
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
``````

@ 2023 Alexandru Ghitza · Created with Zola and PureCSS