function g = CholDiag(u)
% u is n-by-1
% g is an n-by column vector which houses the diagonal of G
% where GG' = I + uu^{T}

n = length(u);
g = zeros(n,1);
% Since g(1)^2 = 1 + u(1)^2...
g(1) = sqrt(1 + u(1)^2);
Denom = 1;
for k=1:n-1
% Since
% [G1 0 ; v' G(k+1,k+1)] [G1 0 ; v' G(k+1,k+1)]' =
% [ I+ww^{T} u(k+1)w ; u(k+1)w' 1+u(k+1)^2]
% where w = u(1:k), it follows that
% G1*G1' = (I + ww')
% G1*v = u(k+1)*w
% G(k+1,k+1)^2 + v'*v = (1+u(k+1)^2)
% and so
% G(k+1,k+1)^2 = (1+u(k+1)^2) - v'*v
% = (1+u(k+1)^2) - (u(k+1)^2)w'*inv(I+ww')w
% By Sherman Morrison,
% inv(I + ww') = I -alfa*ww'
% where alfa = 1/(1+w'*w), so
% w'*inv(I + ww')w = w'*w - (w'w)^2/(1+w'*w) = w'*w/(1+w'*w).
% Thus
% G(k+1,k+1)^2 = 1 + u(k+1)^2(1 - w'*w/(1+w'*w))
% = 1 + u(k+1)^2/(1+w'*w)
%
% Denom = 1 + u(1:k)'*u(1:k)...
Denom = Denom + u(k)^2;
g(k+1) = sqrt(1 + u(k+1)^2/Denom);
end

% Computing w'*w from scratch = -2