16) Linear Algebra

Today¶

Matrices as linear transformations
Polynomial evaluation and fitting
Orthogonality

using Plots
default(linewidth=4, legendfontsize=12)

1. Matrices as linear transformations¶

Linear algebra is the study of linear transformations on vectors, which represent points in a finite dimensional space. The matrix-vector product $y = A x$ is a linear combination of the columns of $A$ . The familiar definition,

y_i = \sum_j A_{i,j} x_j

(1)

can also be viewed as

y = \Bigg[ A_{:,1} \Bigg| A_{:,2} \Bigg| \dotsm \Bigg] \begin{bmatrix} x_1 \\ x_2 \\ \vdots \end{bmatrix} = \Bigg[ A_{:,1} \Bigg] x_1 + \Bigg[ A_{:,2} \Bigg] x_2 + \dotsb .

(2)

That is, a linear combination of the columns of A.

Math and Julia Notation¶

The notation $A_{i,j}$ corresponds to the Julia syntax A[i,j] and the colon : means the entire range (row or column). So $A_{:,j}$ is the $j$ -th column and $A_{i,:}$ is the $i$ -th row.

The corresponding Julia syntax is A[:,j] and A[i,:].

Julia has syntax for row vectors, column vectors, and arrays.

[1 2 3] # a row vector

1×3 Matrix{Int64}:
 1  2  3

[1; 2; 3] # a column vector

3-element Vector{Int64}:
 1
 2
 3

[1. 2 3; 4 5 6] # a 2x3 real matrix
# compare with Python's syntax: np.array([[1, 2, 3], [4, 5, 6]])

2×3 Matrix{Float64}:
 1.0  2.0  3.0
 4.0  5.0  6.0

[1 2; 4 3]

2×2 Matrix{Int64}:
 1  2
 4  3

[1 0; 0 2; 10 3]

3×2 Matrix{Int64}:
  1  0
  0  2
 10  3

[1; 2 + 1im; 3]' # ' is transpose, and for complex-valued matrices is the conjugate transpose ("adjoint")

1×3 adjoint(::Vector{Complex{Int64}}) with eltype Complex{Int64}:
 1+0im  2-1im  3+0im

Implementing multiplication by row¶

function matmult1(A, x)
    m, n = size(A)
    y = zeros(m)
    for i in 1:m # row index first, i.e., for each row
        for j in 1:n # we iterate over the columns
            y[i] += A[i,j] * x[j] # we apply the familiar definition
        end
    end
    y
end

matmult1 (generic function with 1 method)

A = reshape(1.:12, 3, 4) # a 3x4 matrix with the numbers 1:12

3×4 reshape(::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, 3, 4) with eltype Float64:
 1.0  4.0  7.0  10.0
 2.0  5.0  8.0  11.0
 3.0  6.0  9.0  12.0

x = [10., 0, 0, 0]

4-element Vector{Float64}:
 10.0
  0.0
  0.0
  0.0

matmult1(A, x)

3-element Vector{Float64}:
 10.0
 20.0
 30.0

# Dot product
A[2, :]' * x

20.0

function matmult2(A, x)
    m, n = size(A)
    y = zeros(m)
    for i in 1:m # iterate over rows of A
        y[i] = A[i,:]' * x # this way we use the dot product between the transposed row of A and the whole x
    end
    y
end

matmult2(A, x)

3-element Vector{Float64}:
 10.0
 20.0
 30.0

Curiosity: which one is faster? Which one takes more allocations?

@time matmult1(A,x)
@time matmult2(A,x)

  0.000012 seconds (1 allocation: 80 bytes)
  0.000016 seconds (4 allocations: 368 bytes)

3-element Vector{Float64}:
 10.0
 20.0
 30.0

Implementing multiplication by column¶

function matmult3(A, x)
    m, n = size(A)
    y = zeros(m)
    for j in 1:n # iterate over columns of A
        y += A[:, j] * x[j]
    end
    y
end

matmult3(A, x)

3-element Vector{Float64}:
 10.0
 20.0
 30.0

@time matmult3(A, x)

  0.000011 seconds (13 allocations: 1.016 KiB)

3-element Vector{Float64}:
 10.0
 20.0
 30.0

A * x # built-in matrix-vector multiply

3-element Vector{Float64}:
 10.0
 20.0
 30.0

@time A * x  # We'll use this version

  0.000010 seconds (1 allocation: 80 bytes)

3-element Vector{Float64}:
 10.0
 20.0
 30.0

Check the standard operations documentation page.

2. Polynomial evaluation¶

Polynomial evaluation is (continuous) linear algebra¶

We can evaluate polynomials using matrix-vector multiplication. For example,

- 3x + 5x^3 = \Bigg[ 1 \Bigg|\, x \Bigg|\, x^2 \,\Bigg|\, x^3 \Bigg] \begin{bmatrix}0 \\ -3 \\ 0 \\ 5 \end{bmatrix} .

(3)

using Pkg
Pkg.add("Polynomials")

using Polynomials
P(x) = Polynomial(x)

p = [0, -3, 0, 5] # vector of coefficients for the canonical basis
q = [1, 2, 3, 4]
@show f = P(p)
@show g = P(q)
h = f + g
@show h
@show P(p+q) # we can see that the polynomial evaluation as a linear combination (of matrix columns) is a linear transformation!
x = [0., 1, 2]
h.(x) # can be applied to each element of a vector with the dot operator

   Resolving package versions...

  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`

f = P(p) = Polynomial(-3*x + 5*x^3)
g = P(q) = Polynomial(1 + 2*x + 3*x^2 + 4*x^3)
h = Polynomial(1 - x + 3*x^2 + 9*x^3)
P(p + q) =

Polynomial(1 - x + 3*x^2 + 9*x^3)

3-element Vector{Float64}:
  1.0
 12.0
 83.0

plot(h, legend=:bottomright, xlim=(-2, 2))

Loading...

Polynomial evaluation is (discrete) linear algebra¶

V = [one.(x) x x.^2 x.^3] # Vandermonde matrix

3×4 Matrix{Float64}:
 1.0  0.0  0.0  0.0
 1.0  1.0  1.0  1.0
 1.0  2.0  4.0  8.0

V * p + V * q # same as h.(x)

3-element Vector{Float64}:
  1.0
 12.0
 83.0

V * (p + q) # again, matrix multiplication is a linear transformation

3-element Vector{Float64}:
  1.0
 12.0
 83.0

Vandermonde matrices¶

A Vandermonde matrix is one whose columns are polynomials (monomials) evaluated at discrete points.

V(x) = \begin{bmatrix} 1 \Bigg| x \Bigg| x^2 \Bigg| x^3 \Bigg| \dotsb \end{bmatrix}

(4)

function vander(x, k=nothing)
    if isnothing(k)
        k = length(x)
    end
    m = length(x)
    V = ones(m, k)
    for j in 2:k
        V[:, j] = V[:, j-1] .* x
    end
    V
end

vander (generic function with 2 methods)

@show x = LinRange(-1, 1, 50)
V = vander(x, 4)
scatter(x, V, legend=:bottomright, label = ["V_1" "V_2" "V_3" "V_4"])

x = LinRange(-1, 1, 50) = LinRange{Float64}(-1.0, 1.0, 50)

Loading...

Fitting (polynomial interpolation) is linear algebra¶

\underbrace{\begin{bmatrix} 1 \Bigg| x \Bigg| x^2 \Bigg| x^3 \Bigg| \dotsb \end{bmatrix}}_{V(x)} \Big[ p \Big] = \Bigg[ y \Bigg]

(5)

x1 = [-.9, 0.1, .5, .8]
y1 = [1, 2.4, -.2, 1.3]
scatter(x1, y1, markersize=8)

Loading...

V = vander(x1)
@show size(V)
p = V \ y1 # write y1 in the polynomial basis; left-division, read V^{-1} * y1 (like solving p V = y1)
scatter(x1, y1, markersize=8, xlims=(-1, 1))
plot!(P(p), label="P(p)")
plot!(x, vander(x, 4) * p, label="\$ V(x) p\$", linestyle=:dash)

size(V) = (4, 4)

Loading...

Some common terminology¶

The range of $A$ is the space spanned by its columns. This definition coincides with the range of a function $f(x)$ when $f(x) = A x$ .
The (right) nullspace of $A$ is the space of vectors $x$ such that $A x = 0$ .
The rank of $A$ is the dimension of its range.
A matrix has full rank if the nullspace of either $A$ or $A^T$ is empty (only the 0 vector). Equivalently, if all the columns of $A$ (or $A^T$ ) are linearly independent.
A nonsingular (or invertible) matrix is a square matrix of full rank. We call the inverse $A^{-1}$ and it satisfies $A^{-1} A = A A^{-1} = I$ .

Poll 15.1: If $A \in \mathbb{R}^{m\times m}$ , which of these doesn’t belong?¶

$A$ has an inverse, $A^{-1}$
$\operatorname{rank} (A) = m$
$\operatorname{null}(A) = \{0\}$
$A A^T = A^T A$
$\det(A) \ne 0$
$A x = 0$ implies that $x = 0$

Pkg.add("LinearAlgebra")
using LinearAlgebra

A = rand(4, 4)
B = A' * A - A * A' # no. 4! (only valid for orthogonal matrices - more later)
@show B
det(A)

   Resolving package versions...

  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`

B = [-0.7524380040011185 -0.13701483574932205 -0.29457515938763557 -0.06864902392835437; -0.13701483574932205 -0.21478294270764864 -0.23531565443422542 -0.06757256692336011; -0.29457515938763557 -0.23531565443422542 0.3627463157922426 0.5864642755136569; -0.06864902392835437 -0.06757256692336011 0.5864642755136569 0.604474630916525]

0.1701663345605843

What is an inverse?¶

When we write $x = A^{-1} y$ , we mean that $x$ is the unique vector such that $A x = y$ . (It is rare that we explicitly compute a matrix $A^{-1}$ , though it’s not as “bad” as people may have told you.)

A vector $y$ is equivalent to $\sum_i e_i y_i$ where $e_i$ are columns of the identity. Meanwhile, $x = A^{-1} y$ means that we are expressing that same vector $y$ in the basis of the columns of $A$ , i.e., $\sum_i A_{:,i} x_i$ .

using LinearAlgebra
A = rand(4, 4)

4×4 Matrix{Float64}:
 0.601396  0.455792  0.595225   0.770196
 0.482789  0.850813  0.868196   0.267611
 0.912719  0.445778  0.0174431  0.807863
 0.983887  0.303785  0.251255   0.0804293

A \ A # left-division, read A^{-1} * A; notice I (all 1's on the diagonal, and numerically 0's elsewhere)

4×4 Matrix{Float64}:
  1.0           0.0  7.57181e-18  -1.87279e-17
  3.62171e-17   1.0  1.58208e-16   3.41974e-16
 -3.20193e-18  -0.0  1.0          -2.849e-16
  1.42371e-16   0.0  0.0           1.0

inv(A) * A

4×4 Matrix{Float64}:
  1.0           8.68596e-17   8.51652e-17   3.46604e-17
  5.16831e-17   1.0          -4.49128e-17   1.46353e-18
 -5.61439e-18  -5.3878e-17    1.0          -1.28552e-16
  8.08992e-17  -4.39595e-17   8.40137e-17   1.0

3. Orthogonality¶

The inner product

x^T y = \sum_i x_i y_i

(6)

of real vectors (or columns of a matrix) tells us about their magnitude and about the angle. The norm is induced by the inner product,

\lVert x \rVert = \sqrt{x^T x}

(7)

and the angle $\theta$ is defined by

\cos \theta = \frac{x^T y}{\lVert x \rVert \, \lVert y \rVert} .

(8)

Inner products are bilinear, which means that they satisfy some convenient algebraic properties

\begin{split} (x + y)^T z &= x^T z + y^T z \\ x^T (y + z) &= x^T y + x^T z \\ (\alpha x)^T (\beta y) &= \alpha \beta x^T y \\ \end{split} .

(9)

Examples with inner products¶

x = [0, 1]
y = [1, 1]
@show x' * y
@show y' * x;

x' * y = 1
y' * x = 1

ϕ = pi/6
y = [cos(ϕ), sin(ϕ)]
cos_θ = x'*y / (norm(x) * norm(y))
@show cos_θ
@show cos(ϕ-pi/2);

cos_θ = 0.49999999999999994
cos(ϕ - pi / 2) = 0.4999999999999999

Polynomials can be orthogonal too!¶

x = LinRange(-1, 1, 50)
A = vander(x, 4)
M = A * [.5 0 0 0; # 0.5
         0  1 0 0;  # x
         0  0 1 0]' # x^2
# that is, M = [0.5 | x | x^2]
scatter(x, M, label = ["M_1" "M_2" "M_3"])
plot!(x, 0*x, label=:none, color=:black)

Loading...

Which inner product will be zero?
- Which functions are even and odd?

Polynomial inner products¶

M[:,1]' * M[:,2]

-2.220446049250313e-16

M[:,1]' * M[:,3]

8.673469387755102

M[:,2]' * M[:,3]

-4.440892098500626e-16

Today¶

1. Matrices as linear transformations¶

Math and Julia Notation¶

Implementing multiplication by row¶

Implementing multiplication by column¶

2. Polynomial evaluation¶

Polynomial evaluation is (continuous) linear algebra¶

Polynomial evaluation is (discrete) linear algebra¶

Vandermonde matrices¶

Fitting (polynomial interpolation) is linear algebra¶

Some common terminology¶

Poll 15.1: If A∈Rm×mA \in \mathbb{R}^{m\times m}A∈Rm×m, which of these doesn’t belong?¶

What is an inverse?¶

3. Orthogonality¶

Examples with inner products¶

Polynomials can be orthogonal too!¶

Polynomial inner products¶

Poll 15.1: If $A \in \mathbb{R}^{m\times m}$ , which of these doesn’t belong?¶