##
using LinearAlgebra
using Random

##
function least_squares_eliminate(A,b)
  a = A[:,1]
  na = norm(a)
  a = a/na
  if size(A,2) == 1
    return [a'*b]/na
  end
  y = least_squares_eliminate(A[:,2:end]-a*a'*A[:,2:end], b - a*a'*b)
  γ = a'*(b - A[:,2:end]*y)/na
  return pushfirst!(y,γ)
end
Random.seed!(2)
A = randn(5,3)
b = randn(5)
xtrue = A\b
display([least_squares_eliminate(A,b), xtrue])
norm(xtrue-least_squares_eliminate(A,b))

2-element Array{Array{Float64,1},1}:
 [-0.5294972693999612, -0.006905042606022434, 0.15718401182671496]
 [-0.5294972693999611, -0.006905042606022456, 0.15718401182671493]

1.1663999956178395e-16

## Let's see if my soln is the same
function least_squares_transform(A)
  T = Matrix(1.0I,size(A,1),size(A,1))
  Z = Matrix(1.0I,size(A,1),size(A,1))
  R = Matrix(1.0I,size(A,1),size(A,1))
  Q = Matrix(0.0I,size(A,1),size(A,2))
  for i=1:size(A,2)
    a = A[:,i]
    na = norm(a)
    a = a/na
    Q[:,i] = a
    R = R*(Matrix(1.0I,size(A,1),size(A,1)) - a*a')
    T = (Matrix(1.0I,size(A,1),size(A,1)) - a*a')*T
    A = (Matrix(1.0I,size(A,1),size(A,1)) - a*a')*A
    Z = (Matrix(1.0I,size(A,1),size(A,1)) - 2*a*a')*Z
  end
  return T,R,Q,Z
end
T,R,Q,Z = least_squares_transform(A)
display(T*b)

5-element Array{Float64,1}:
  0.15001460522896287 
  0.09347796911404205 
 -0.6802657944831366  
  0.25281562752716547 
  0.028994342988093974