using LinearAlgebra

# 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

# 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₄]

# 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₄]

#check that Q has orthonormal columns
Q'Q

Q'Q - I

Q'Q ≈ I

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

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

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

R = F.R

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

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

R # Recognize this matrix?

Q2*R ≈ A

b = rand(6)

A \ b

inv(A'A) * A'b

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

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