Getting started

https://julialang.org/

• Recommended: juliaup

Getting started

https://www.cs.cornell.edu/~bindel/nmds/

• See chapter 2 (Julia programming)
• You will not need all of the macro stuff!

Jupyter, Pluto, Quarto

For noodling around, notebook environments are great

• Jupyter
  • Many languages: Julia, Python, R, and others
  • Google Colab is built on Jupyter
• Pluto
  • Julia only
  • Plain text format, notebooks are Julia files
  • Reactive notebooks (vs “hidden state”)
  • Will use this for assignments and such

For presentation, I use Quarto

Starting Pluto

After installing Pluto:

using Pluto
Pluto.run()

For most of this class

Start with

using LinearAlgebra
using Plots

Command line demo!

Linear algebra

Representation of

• A linear map \(\mathbb{R}^n \rightarrow \mathbb{R}^m\)
• More general linear maps in some basis
• Maybe other things! Ex: quadratic forms
  • More on this next week

Matrix algebra

\[A \in \mathbb{R}^{m \times n}\]

• An array of real numbers
• Two indices: \(a_{ij}\) is row \(i\), column \(j\)
• An element of a vector space \(\mathbb{R}^{m \times n}\)

Perspectives

\[A \in \mathbb{R}^{m \times n}\]

\[A = \begin{bmatrix} a_{1,:} \\ \vdots \\ a_{m,:} \end{bmatrix}\]

\[A = \begin{bmatrix} a_{:,1} & \ldots & a_{:,n} \end{bmatrix}\]

Julia

A = [1.0 2.0; 3.0 5.0]
B = [7.0 8.0; 9.0 12.0]
x = [1.0, 2.0]
y = [3.0, 4.0]

Commandment 1

Matrix \(\times\) vector \(=\) linear combination of matrix columns.

• Think \(A = \begin{bmatrix} a_1 & \ldots & a_n \end{bmatrix}\)
• Then \(Ax = \sum_{j=1}^n a_j x_j\)

function matvec_col(A, x)
    m, n = size(A)
    y = zeros(m)
    for j = 1:n
        y += A[:,j]*x[j]
    end
    y
end

matvec_col(A, x) ≈ A*x

true

Commandment 2

Inner product \(=\) sum of products of corresponding elements.

• This is for the standard inner product for a real vector space: \[x \cdot y = \sum_{j=1}^n x_j y_j\]

d1 = sum(xi*yi for (xi,yi) in zip(x,y))
d2 = dot(x, y)
d3 = y'*x

d1 ≈ d3, d2 ≈ d3

(true, true)

Commandment 3

Consider \(u, v, w \in \mathbb{R}^n\) where \(n = 10^6\), and let \(y = Aw\) for \(A = uv^T\)

y = u*v'*w   # Slow!

A = u*v'
y = A*w      # Also slow!

y = u*(v'*w) # Much better

Commandment 4

Matrix \(\times\) diagonal \(=\) scaling of the matrix columns.

A * Diagonal(x), A .* x'

([1.0 4.0; 3.0 10.0], [1.0 4.0; 3.0 10.0])

Commandment 5

Diagonal \(\times\) matrix \(=\) scaling of the matrix rows.

Diagonal(x) * A, x .* A

([1.0 2.0; 6.0 10.0], [1.0 2.0; 6.0 10.0])

Commandment 6

Never form an explicit diagonal matrix.

Diagonal(x), diagm(x)

(Diagonal([1.0, 2.0]), [1.0 0.0; 0.0 2.0])

Commandment 7

Never form an explicit rank one matrix.

• Storage for \(u, v \in \mathbb{R}^n\): \(O(n)\)
• Computing \(u (v^T w)\): \(O(n)\) flops
• Forming \(uv^T\): \(O(n^2)\) flops and storage

Commandment 8

Matrix \(\times\) matrix \(=\) collection of matrix-vector products.

• If \(B = \begin{bmatrix} b_1 & \ldots & b_n \end{bmatrix}\) then \(AB = \begin{bmatrix} Ab_1 & \ldots & Ab_n \end{bmatrix}\)

function matmul_col(A, B)
    m, p = size(A)
    p, n = size(B)
    C = zeros(m, n)
    for j = 1:n
        C[:,j] = A*B[:,j]
    end
    C
end

matmul_col(A, B) ≈ A*B

true

Commandment 9

Matrix \(\times\) matrix \(=\) collection of dot products.

function matmul_dots(A, B)
    m, p = size(A)
    p, n = size(B)
    C = zeros(m, n)
    for i = 1:m
        for j = 1:n
            C[i,j] = A[i,:]'*B[:,j]
        end
    end
    C
end

matmul_dots(A, B) ≈ A*B

true

Commandment 10

Matrix \(\times\) matrix \(=\) sum of rank one matrices.

function matmul_outer(A, B)
    m, p = size(A)
    p, n = size(B)
    C = zeros(m, n)
    for k = 1:p
        C += A[:,k]*B[k,:]'
    end
    C
end

matmul_outer(A, B) ≈ A*B

true

Commandment 11

Matrix \(\times\) matrix \(\implies\) linearly combine rows from the second matrix.

function matmul_rows(A, B)
    m, p = size(A)
    p, n = size(B)
    C = zeros(m, n)
    for i = 1:m
        C[i,:] = A[i,:]'*B
    end
    C
end

matmul_rows(A, B) ≈ A*B

true

Commandment 12

Matrix \(\times\) matrix \(\implies\) linearly combine columns from the first matrix.

This is really the same as saying

\[A \begin{bmatrix} b_1 & \ldots & b_n \end{bmatrix} = \begin{bmatrix} Ab_1 & \ldots & Ab_n \end{bmatrix}\]

Example

\[M = \begin{bmatrix} A & b \\ c^T & d \end{bmatrix} = \begin{bmatrix} a_{11} & a_{12} & b_1 \\ a_{21} & a_{22} & b_2 \\ c_1 & c_2 & d \end{bmatrix}\]

Julia

b = [7; 11]
c = [13; 14]
d = 42
M = [A b; c' d]

3×3 Matrix{Float64}:
  1.0   2.0   7.0
  3.0   5.0  11.0
 13.0  14.0  42.0

Think blocks

M * [x; 10] ≈
[A*x + b*10; c'*x + d*10]

true

What about

[I b; c' d] ≈
[Diagonal(x) b; c' d]

false

Dense square matrix

Arand = rand(5,5)

5×5 Matrix{Float64}:
 0.422307   0.877403   0.454838  0.611675  0.936974
 0.397803   0.676438   0.704521  0.315096  0.92584
 0.304245   0.0519322  0.452553  0.152442  0.0138345
 0.690428   0.212928   0.377369  0.302685  0.891633
 0.0312779  0.803481   0.951678  0.77644   0.772113

Upper triangular

triu(Arand)

5×5 Matrix{Float64}:
 0.422307  0.877403  0.454838  0.611675  0.936974
 0.0       0.676438  0.704521  0.315096  0.92584
 0.0       0.0       0.452553  0.152442  0.0138345
 0.0       0.0       0.0       0.302685  0.891633
 0.0       0.0       0.0       0.0       0.772113

Strict upper triangle

triu(Arand,1)

5×5 Matrix{Float64}:
 0.0  0.877403  0.454838  0.611675  0.936974
 0.0  0.0       0.704521  0.315096  0.92584
 0.0  0.0       0.0       0.152442  0.0138345
 0.0  0.0       0.0       0.0       0.891633
 0.0  0.0       0.0       0.0       0.0

Upper triangle plus a bit

triu(Arand,1)

5×5 Matrix{Float64}:
 0.0  0.877403  0.454838  0.611675  0.936974
 0.0  0.0       0.704521  0.315096  0.92584
 0.0  0.0       0.0       0.152442  0.0138345
 0.0  0.0       0.0       0.0       0.891633
 0.0  0.0       0.0       0.0       0.0

Upper triangular

UpperTriangular(Arand)

5×5 UpperTriangular{Float64, Matrix{Float64}}:
 0.422307  0.877403  0.454838  0.611675  0.936974
  ⋅        0.676438  0.704521  0.315096  0.92584
  ⋅         ⋅        0.452553  0.152442  0.0138345
  ⋅         ⋅         ⋅        0.302685  0.891633
  ⋅         ⋅         ⋅         ⋅        0.772113

Upper triangular

tril(Arand)

5×5 Matrix{Float64}:
 0.422307   0.0        0.0       0.0       0.0
 0.397803   0.676438   0.0       0.0       0.0
 0.304245   0.0519322  0.452553  0.0       0.0
 0.690428   0.212928   0.377369  0.302685  0.0
 0.0312779  0.803481   0.951678  0.77644   0.772113

Upper triangular

LowerTriangular(Arand)

5×5 LowerTriangular{Float64, Matrix{Float64}}:
 0.422307    ⋅          ⋅         ⋅         ⋅ 
 0.397803   0.676438    ⋅         ⋅         ⋅ 
 0.304245   0.0519322  0.452553   ⋅         ⋅ 
 0.690428   0.212928   0.377369  0.302685   ⋅ 
 0.0312779  0.803481   0.951678  0.77644   0.772113

Spy plots

spy(triu(rand(100,100)))

Sparse matrix

For \(A \in \mathbb{R}^{m \times n}\)

• Define \(\operatorname{nnz}(A) =\) number of nonzero entries
• Say \(A\) is sparse if \(\operatorname{nnz}(A) \ll mn\)
• Example: diagonal matrices

Sparse matrix

Why is \(A\) sparse interesting?

• Maybe more compact storage (sparse formats)
• Maybe faster operations (e.g. matrix-vector multiply is \(O(\operatorname{nnz}(A))\))

Will discuss these points more later on

Data-sparse matrix

A dense matrix can still have compact storage / fast multiply!

Examples:

• Rank one matrix \(A = uv^T\)
• Upper triangle of a rank one matrix
• Topelitz matrix (constant diagonal entries)