Robot Revolution

What is "Linear Algebra"?

There are plenty of technical definitions, on [Wikipedia][linalg-wiki], in textbooks and on excellent websites such as [3Blue1Brown][3blue1brown], but for our purposes we can be more informal.

The Julia (and Python) maintainers love this sort of thing. Exercism management, not so much.

This Concept will be restricted to the simpler aspects of a huge subject, but be warned: this is, inevitably, quite a mathematical concept.

What is a Vector?

It depends who you ask, and how you want to visualize something rather abstract.

• In Julia, Vector{T} is just an alias for Array{T, 1}, and the eltype T can be anything.
• In (much of) Physics, a vector is an arrow with a length and direction in N-dimensional space, but no fixed position.
• In Linear Algebra, the physicists' arrow is fixed in space, with its tail at the origin. This makes the shaft of the arrow redundant, so we can represent the vector as a point in N-dimensional space, with the N elements representing distance from the origin along each of the N axes.

For present purposes, ignore vectors of strings or characters. We will be working with numerical types in this Concept: Int, Float, or Complex.

What is a matrix?

Again, you will find a variety of answers.

• A rectangular 2-D array of numbers (with no ragged edges). A square matrix is a common special case.
• A linear combination of column vectors, stacked side by side.
• A transformation that can be applied to a vector, just as a function is applied to other input types.

Some types of square matrix are common enough to have special names.

Diagonal matrix: All non-zero entries are on the main diagonal (top-left to bottom-right).

julia> [1 0 0; 0 2 0; 0 0 3]
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

# a convenient shortcut
julia> diagm(1:3)
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

Identity matrix: A diagonal matrix with only ones on the diagonal (for reasons that will become clearer later). Often abbreviated as I.

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

Upper triangular: Non-zero values on and above the diagonal, zeros below. Lower triangular you can guess.

Transposing a matrix

If you really want to swap rows with columns, the function permutedims() will do this and give you a new matrix.

This involves copying, which is slow and memory-hungry when working with large matrices.

For Linear Algebra purposes, the transpose() function is more useful, as it quickly creates a lazy wrapper around the original matrix.

Even more useful is the adjoint() function, which also flips the sign of the imaginary part in any complex numbers. This operation is common enough that we can just add an apostrophe ' to the variable name to create the adjoint.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

# new, full copy
julia> permutedims(m)
3×2 Matrix{Int64}:
 1  4
 2  5
 3  6

# lazy version
julia> transpose(m)
3×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6

julia> mc = [1+2im 2+3im; 3+2im 1+2im]
2×2 Matrix{Complex{Int64}}:
 1+2im  2+3im
 3+2im  1+2im

# lazy conjugate transpose
julia> adjoint(mc)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

# syntactic sugar for adjoint
julia> mc'
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

Multiplication

The rest of this document includes functionality from the LinearAlgebra module. A using LinearAlgebra line will bring it into the namespace, but we will not keep repeating this in the examples (too much visual clutter).

Vectors, element-wise product

This was discussed in the [Vector Operations][vector-ops] Concept. The operator is .*, which operates on the input vectors pair-wise, giving an output the same size and type as the inputs.

julia> [1, 2] .* [3, 4]
2-element Vector{Int64}:
 3
 8

Vectors, dot product

This extremely common operation is written in textbooks as u ⋅ v, and called the "dot" product.

A dot product is equivalent to the sum of the element-wise product.

Two vectors can be multiplied with the usual * operator, but only if the left vector is an adjoint: conveniently written u' * v. This should become clearer in the later section on matrix multiplication.

The raised (center) dot is available in Julia (entered \cdot tab) as syntactic sugar for the dot() function. There is no need to specify the adjoint, as this detail is handled automatically.

julia> using LinearAlgebra

# with element-wise syntax
julia> sum([1, 2] .* [3, 4])
11

# with u' * v syntax
julia> [1, 2]' * [3, 4]
11

# with dot()
julia> dot([1, 2], [3, 4])
11

# with \cdot syntax
julia> [1, 2] ⋅ [3, 4]
11

For complex-valued vectors, the left vector needs to be the conjugate (sign flipped on the imaginary part). The dot function does this automatically, the u' syntax does it explicitly, but sum(u .* v) will fail if u and v are complex.

Norms

How "big" is a vector?

The norm is an attempt to capture this, by reducing the vector to an appropriate scalar.

A whole family of norms are defined, but by far the most common is the 2-norm, which is √(v ⋅ v).

This root-mean-square operation is the Pythagorean distance from the origin (in N-dimensional space). If we visualize the vector as an arrow, with its tail at the origin, the 2-norm is the length of the arrow.

Any p-norm can be calculated by providing p as the 2nd argument. The 1-norm is sometimes useful: it is simply the sum of the absolute values of the entries (so very quick and easy to calculate).

# defaults to the 2-norm
julia> norm([1, 2, 3])
3.7416573867739413
# the 1-norm
julia> norm([1, -2, 3], 1)
6.0

Matrix multiplication

It sometimes seems as though every applied mathematician in the world has spent much of the last 80 years turning every calculation into a series of matrix multiplications.

This is an operation that computers are very good at:

• It is highly repetitive.
• It can be parallelized efficiently.
• A lot of specialized hardware has been developed to make it faster, including the GPU that is probably integrated into your laptop.

The details are fairly simple, though not very intuitive at first sight.

Consider a matrix A multiplying a vector v to give an output w (by Linear algebra convention, we will uses uppercase letters for matrices, lowercase for vectors).

The top row of A is dotted with v to get the first element of w, the second row gives the second element, and so on down.

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4
julia> v = [5, 6]
2-element Vector{Int64}:
 5
 6
julia> A * v
2-element Vector{Int64}:
 17  # equals [1, 2] ⋅ [5, 6]
 39  # equals [3, 4] ⋅ [5, 6]

A matrix * matrix multiplication extends this across the columns of the matrix on the right.

For C = A * B, we can consider that A multiplies each column of B to give the corresponding column of C: a series of matrix-vector multiplications.

Equivalently, we could say that the first row of A is dotted with each column of B to give the top row of C, the second row gives the second row, and so on down.

julia> A
2×2 Matrix{Int64}:
 1  2
 3  4

julia> B = [5 7; 6 8]
2×2 Matrix{Int64}:
 5  7
 6  8

# left column is the same as A*v previously
julia> A * B
2×2 Matrix{Int64}:
 17  23
 39  53
 
julia> B * A
2×2 Matrix{Int64}:
 26  38
 30  44

As shown in the example above, matrix multiplication does not commute: there is no simple relationship between A*B and B*A.

Matrix multiplication is probably hard to visualize for anyone new to it, just by reading the words. YouTube has lots of videos demonstrating it graphically, so search for "matrix multiplication" and choose one with your preferred style, detail level and language.

Dimensions

The dot product of two vectors relies on them having equal length.

By extension, for matrix multiplication the number of columns of the left matrix must match the number of rows of the right matrix.

Expressing the sizes as (nrows, ncols) tuples, as output by size(A), we have (a, b) * (b, c) -> (a, c). The "inner" dimensions, here b and b, are compatible with taking dot products. The "outer" dimensions, here a and c, determine the dimensions of the output.

An example with rectangular matrices:

julia> D = reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> E = reshape(1:12, 3, 4)
3×4 reshape(::UnitRange{Int64}, 3, 4) with eltype Int64:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> D * E
2×4 Matrix{Int64}:
 22  49   76  103
 28  64  100  136

# (3, 4) * (2, 3) not possible
julia> E * D
ERROR: DimensionMismatch: matrix A has axes (Base.OneTo(3),Base.OneTo(4)), matrix B has axes (Base.OneTo(2),Base.OneTo(3))

Multiplying pairs of vectors, matrix-style, has two possibilities.

Conventionally, we would use u' * v as an equivalent to the dot product. The dimensions are (1, 3) * (3, 1), and Julia simplifies the (1, 1) output to a scalar (unlike, for example, R).

Alternatively, we could use u * v', with dimensions (3, 1) * (1, 3) -> (3, 3), to multiply all possible pairs of elements and expand the vectors to a matrix.

julia> u = [1, 2]
2-element Vector{Int64}:
 1
 2
julia> v = [3, 4]
2-element Vector{Int64}:
 3
 4
julia> u' * v
11
julia> u * v'
2×2 Matrix{Int64}:
 3  4
 6  8

u' * v is sometimes called the inner product.

Rotations

One particularly common type of matrix multiplication involves rotation matrices.

In 2D, there is a relatively simple matrix to rotate a vector anti-clockwise by θ radians.

julia> rot2d(θ, vec) = [cos(θ) -sin(θ); sin(θ) cos(θ)] * vec
rot2d (generic function with 1 method)

# unit vector in the x direction
julia> i_hat = [1, 0]
2-element Vector{Int64}:
 1
 0

# rotate 45 degrees
julia> rot2d(π/4, i_hat)
2-element Vector{Float64}:
 0.7071067811865476
 0.7071067811865475

# rotate 90 degrees -> unit vector in the y direction, j_hat
julia> rot2d(π/2, i_hat)
2-element Vector{Float64}:
 6.123233995736766e-17  # zero, within numerical error
 1.0