# 2022-02-11 Gram-Schmidt

## Contents

# 2022-02-11 Gram-Schmidt¶

## Last time¶

Projections

Rotations

Reflections

## Today¶

Revisit projections, rotations, and reflections

Constructing orthogonal bases

QR factorization

```
using LinearAlgebra
using Plots
using Polynomials
default(linewidth=4, legendfontsize=12)
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
```

```
┌ Info: Precompiling Polynomials [f27b6e38-b328-58d1-80ce-0feddd5e7a45]
└ @ Base loading.jl:1423
```

```
vander (generic function with 2 methods)
```

# 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
scatter(x, M)
plot!(x, 0*x, label=:none, color=:black)
```

Which inner product will be zero?

Which functions are even and odd?

## Pairwise inner product¶

The pairwise inner products between two sets of vectors can be expressed by collecting the sets as columns in matrices and writing \(A = X^T Y\) where \(A_{i,j} = x_i^T y_j\). It follows from this definition that

```
M' * M
```

```
3×3 Matrix{Float64}:
12.5 -1.11022e-15 8.67347
-1.11022e-15 17.3469 -2.22045e-16
8.67347 -2.22045e-16 10.8272
```

# Normalization and orthogonalization¶

```
q1 = M[:,1]
q1 /= norm(q1) # normalize q1
Q = [q1 M[:, 2:end]]
scatter(x, Q)
```

```
Q' * Q
```

```
3×3 Matrix{Float64}:
1.0 -1.94289e-16 2.45323
-1.94289e-16 17.3469 -2.22045e-16
2.45323 -2.22045e-16 10.8272
```

# Orthogonal matrices¶

If \(x^T y = 0\) then we say \(x\) and \(y\) are **orthogonal** (or “\(x\) is orthogonal to \(y\)”).
A vector is said to be **normalized** if \(\lVert x \rVert = 1\).
If \(x\) is orthogonal to \(y\) and \(\lVert x \rVert = \lVert y \rVert = 1\) then we say \(x\) and \(y\) are **orthonormal**.
A square matrix with orthonormal columns is said to be an **orthogonal matrix**.
We typically use \(Q\) or \(U\) and \(V\) for matrices that are known/constructed to be orthogonal.
Orthogonal matrices are always full rank – the columns are linearly independent.
The inverse of a orthogonal matrix is its transpose:

# Gram-Schmidt orthogonalization¶

Suppose we’re given some vectors and want to find an orthogonal basis for their span.

# A naive algorithm¶

```
function gram_schmidt_naive(A)
m, n = size(A)
Q = zeros(m, n)
R = zeros(n, n)
for j in 1:n
v = A[:,j]
for k in 1:j-1
r = Q[:,k]' * v
v -= Q[:,k] * r
R[k,j] = r
end
R[j,j] = norm(v)
Q[:,j] = v / R[j,j]
end
Q, R
end
```

```
gram_schmidt_naive (generic function with 1 method)
```

```
A = vander(x, 4)
Q, R = gram_schmidt_naive(A)
@show norm(Q' * Q - I)
@show norm(Q * R - A);
```

```
norm(Q' * Q - I) = 6.609184312080184e-16
norm(Q * R - A) = 2.500206551829073e-16
```

# What do orthogonal polynomials look like?¶

```
A = vander(x, 5)
Q, R = gram_schmidt_naive(A)
plot(x, Q, marker=:auto)
```

What happens if we use more than 50 values of \(x\)? Is there a continuous limit?

# Theorem¶

## Every full-rank \(m\times n\) matrix (\(m \ge n\)) has a unique reduced \(Q R\) factorization with \(R_{j,j} > 0\)¶

The algorithm we’re using generates this matrix due to the line:

```
R[j,j] = norm(v)
```

# Solving equations using \(QR = A\)¶

If \(A x = b\) then \(Rx = Q^T b\).

```
x1 = [-0.9, 0.1, 0.5, 0.8] # points where we know values
y1 = [1, 2.4, -0.2, 1.3]
scatter(x1, y1)
Q, R = gram_schmidt_naive(vander(x1, 3))
p = R \ (Q' * y1)
plot!(x, vander(x, 3) * p)
```

# How accurate is it?¶

```
m = 20
x = LinRange(-1, 1, m)
A = vander(x, m)
Q, R = gram_schmidt_naive(A)
@show norm(Q' * Q - I)
@show norm(Q * R - A)
```

```
norm(Q' * Q - I) = 1.073721107832196e-8
norm(Q * R - A) = 8.268821431611631e-16
```

```
8.268821431611631e-16
```

# A variant with more parallelism¶

```
function gram_schmidt_classical(A)
m, n = size(A)
Q = zeros(m, n)
R = zeros(n, n)
for j in 1:n
v = A[:,j]
R[1:j-1,j] = Q[:,1:j-1]' * v
v -= Q[:,1:j-1] * R[1:j-1,j]
R[j,j] = norm(v)
Q[:,j] = v / norm(v)
end
Q, R
end
```

```
gram_schmidt_classical (generic function with 1 method)
```

```
m = 10
x = LinRange(-1, 1, m)
A = vander(x, m)
Q, R = gram_schmidt_classical(A)
@show norm(Q' * Q - I)
@show norm(Q * R - A)
```

```
norm(Q' * Q - I) = 6.339875256299394e-11
norm(Q * R - A) = 1.217027619812654e-16
```

```
1.217027619812654e-16
```