Gram–Schmidt orthogonalization¶

Chapter 4.4 illustrates a hand technique for computing orthonormal vectors q₁,q₂,… from arbitrary vectors a,b,… with the property that the first k vectors in the original set span the same subspace as the orthonormal set, and this is true for k=1,2,3,...

We will move this hand technique to the computer in this notebook. Some of you will notice that on the computer one can combine operations in a simpler block fashion.

In [1]:

using LinearAlgebra

In [2]:

# start with four arbitrary independent vectors in ℝᵐ
# with random entries from 1 to 10.
m = 6
a₁ = rand(1:10,m)
a₂ = rand(1:10,m)
a₃ = rand(1:10,m)
a₄ = rand(1:10,m)
A = [a₁ a₂ a₃ a₄] # show them as the columns of a 6×4 matrix A

Out[2]:

6×4 Matrix{Int64}:
  5  10  6  8
 10   6  5  1
  2   1  7  3
  7   5  4  8
  1   3  8  3
  1   8  9  4

In [3]:

# The vₖ are vectors, but they are all orthogonal and
#span([v₁]) = span([a₁])
#span([v₁ v₂]) = span([a₁ a₂])
#span([v₁ v₂ v₃]) = span([a₁ a₂ a₃] )
#span([v₁ v₂ v₃ v₄]) = span([a₁ a₂ a₃ a₄])
v₁ = a₁
v₂ = a₂ - v₁*(v₁'a₂)/(v₁'v₁)
v₃ = a₃ - v₁*(v₁'a₃)/(v₁'v₁) - v₂*(v₂'a₃)/(v₂'v₂)
v₄ = a₄ - v₁*(v₁'a₄)/(v₁'v₁) - v₂*(v₂'a₄)/(v₂'v₂) - v₃*(v₃'a₄)/(v₃'v₃)

# gather into a matrix V with orthogonal but *not* orthonormal columns
V = [v₁ v₂ v₃ v₄]

Out[3]:

6×4 Matrix{Float64}:
  5.0   5.61111   -3.16215     1.37162
 10.0  -2.77778   -0.0979465  -3.70534
  2.0  -0.755556   6.16936     1.22464
  7.0  -1.14444   -0.324354    4.20193
  1.0   2.12222    5.22283     0.0756904
  1.0   7.12222    1.49913    -1.74319

In [4]:

# now we normalize
q₁ = normalize(v₁)
q₂ = normalize(v₂)
q₃ = normalize(v₃)
q₄ = normalize(v₄);

# Gather into a matrix Q with orthonormal columns
Q = [q₁ q₂ q₃ q₄]

Out[4]:

6×4 Matrix{Float64}:
 0.372678    0.571756   -0.358733    0.223061
 0.745356   -0.283047   -0.0111116  -0.602583
 0.149071   -0.0769889   0.699888    0.199157
 0.521749   -0.116616   -0.0367966   0.683342
 0.0745356   0.216248    0.592508    0.0123092
 0.0745356   0.725734    0.170071   -0.283488

In [5]:

#check that Q has orthonormal columns
Q'Q

Out[5]:

4×4 Matrix{Float64}:
  1.0          -9.21828e-17   8.99478e-17   1.35877e-16
 -9.21828e-17   1.0           1.0142e-17    1.42552e-16
  8.99478e-17   1.0142e-17    1.0          -2.09229e-16
  1.35877e-16   1.42552e-16  -2.09229e-16   1.0

In [6]:

Q'Q - I

Out[6]:

4×4 Matrix{Float64}:
  0.0          -9.21828e-17   8.99478e-17   1.35877e-16
 -9.21828e-17   0.0           1.0142e-17    1.42552e-16
  8.99478e-17   1.0142e-17    0.0          -2.09229e-16
  1.35877e-16   1.42552e-16  -2.09229e-16   2.22045e-16

In [7]:

Q'Q ≈ I

Out[7]:

true

In [8]:

# compare to what happens if we didn't normalize:
V'V # = diagonal matrix (orthogonal columns, but not orthonormal)

Out[8]:

4×4 Matrix{Float64}:
 180.0          -1.42109e-14   1.33227e-14   1.28786e-14
  -1.42109e-14  96.3111        9.32625e-15   6.94211e-15
   1.33227e-14   9.32625e-15  77.7003       -1.16249e-14
   1.28786e-14   6.94211e-15  -1.16249e-14  37.8113

In [9]:

# What does this triangular structure say?
round.(Q'A, digits=5)

Out[9]:

4×4 Matrix{Float64}:
 13.4164  11.7766   10.3605   8.86974
 -0.0      9.81382   9.2715   6.67879
  0.0      0.0       8.81478  1.38213
  0.0      0.0       0.0      6.14909

QR factorization¶

How do we do all this at once on a computer? We ask the computer to factor the matrix as $QR$ (orthonormal columns times upper triangular).

In [10]:

F = qr(A)  # returns a "factorization object" that stores both Q (implicitly) and R

Out[10]:

LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}}
Q factor:
6×6 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}:
 -0.372678    0.571756    0.358733   -0.223061   -0.0840586  -0.590504
 -0.745356   -0.283047    0.0111116   0.602583    0.0189054  -0.0272178
 -0.149071   -0.0769889  -0.699888   -0.199157   -0.619405   -0.242243
 -0.521749   -0.116616    0.0367966  -0.683342    0.131189    0.478182
 -0.0745356   0.216248   -0.592508   -0.0123092   0.74464    -0.204877
 -0.0745356   0.725734   -0.170071    0.283488   -0.192913    0.566789
R factor:
4×4 Matrix{Float64}:
 -13.4164  -11.7766   -10.3604   -8.86974
   0.0       9.81382    9.2715    6.67879
   0.0       0.0       -8.81478  -1.38213
   0.0       0.0        0.0      -6.14909

In [11]:

R = F.R

Out[11]:

4×4 Matrix{Float64}:
 -13.4164  -11.7766   -10.3604   -8.86974
   0.0       9.81382    9.2715    6.67879
   0.0       0.0       -8.81478  -1.38213
   0.0       0.0        0.0      -6.14909

In [12]:

Q2 = Matrix(F.Q)  # extract the "thin" QR factor you would get from Gram–Schmidt

Out[12]:

6×4 Matrix{Float64}:
 -0.372678    0.571756    0.358733   -0.223061
 -0.745356   -0.283047    0.0111116   0.602583
 -0.149071   -0.0769889  -0.699888   -0.199157
 -0.521749   -0.116616    0.0367966  -0.683342
 -0.0745356   0.216248   -0.592508   -0.0123092
 -0.0745356   0.725734   -0.170071    0.283488

In [13]:

round.(Q'Q2, digits=5) # almost I, up to signs

Out[13]:

4×4 Matrix{Float64}:
 -1.0  -0.0  -0.0   0.0
  0.0   1.0   0.0  -0.0
 -0.0   0.0  -1.0   0.0
 -0.0  -0.0  -0.0  -1.0

In [14]:

R # Recognize this matrix?

Out[14]:

4×4 Matrix{Float64}:
 -13.4164  -11.7766   -10.3604   -8.86974
   0.0       9.81382    9.2715    6.67879
   0.0       0.0       -8.81478  -1.38213
   0.0       0.0        0.0      -6.14909

In [15]:

Q2*R ≈ A

Out[15]:

true

In [16]:

b = rand(6)

Out[16]:

6-element Vector{Float64}:
 0.5450097629781777
 0.8599801580391012
 0.7036387925825908
 0.7553540639899048
 0.7234080262946185
 0.14725528162868073

In [17]:

A \ b

Out[17]:

4-element Vector{Float64}:
  0.08995466028300682
 -0.0768023071730735
  0.07513430689992019
  0.03688677949256471

In [18]:

inv(A'A) * A'b

Out[18]:

4-element Vector{Float64}:
  0.08995466028300675
 -0.07680230717307351
  0.07513430689992023
  0.03688677949256475

In [19]:

R \ (Q2'b)[1:4]

Out[19]:

4-element Vector{Float64}:
  0.0899546602830068
 -0.07680230717307351
  0.07513430689992022
  0.03688677949256465

In [20]:

F \ b  # the factorization object F can be used directly for a least-square solve

Out[20]:

4-element Vector{Float64}:
  0.08995466028300678
 -0.07680230717307349
  0.07513430689992023
  0.03688677949256459

In [ ]: