function sn = SigmaMinEst(U)
% U is an n-by-n nonsingular upper triangular matrix.
% sn is an estimate of its smallest singular value

% Idea. Construct unit n-vector d so that the solution
% to Uy = d is large in norm. Then set sn = 1/norm(y)

[n,n] = size(U);
d = zeros(n,1); d(n) = 1;
y = zeros(n,1); y(n) = 1/U(n,n);
gamma = abs(y(n));
for k=n-1:-1:1
% Determine s = sin(theta) and c = cos(theta) so that if w(k:n) solves
% U(k:n,k:n)w(k:n) = [c;s*d(k+1:n)]
% then norm(w(k:n)) is maximized. Note that since
% U(k+1:n,k+1)*y(k+1:n) = d(k+1:n)
% then
% w(k+1:n) = U(k+1:n,k+1:n)\(s*d(k+1:n)) = s*y(k+1:n)
% Thus,
% w(k) = (c - U(k,k+1:n)*w(k+1:n))/U(k,k)
% = (c - s*U(k,k+1:n)*y(k+1:n))/U(k,k)

% It follows that
% norm(w(k:n))^2 =((c-s*U(k,k+1:n)*y(k+1:n))/U(k,k))^2 + s^2norm(y(k+1:n))^2
% = norm(M*[c;s])^2
% where
alfa = 1/U(k,k); beta = -U(k,k+1:n)*y(k+1:n)/U(k,k);
M = [ alfa beta ; 0 gamma];
% The SVD of M can be used to find the optimal c and s...
[U1,S,V] = svd(M);
c = V(1,1); s = V(2,1);
% Update y and its norm...
y(k:n) = [alfa*c + beta*s; s*y(k+1:n)];
gamma = sqrt(y(k)^2 + (s*gamma)^2);
end
sn = 1/gamma;