A = rand(3,3)

eigA = eigvals(A)

# 1/eigenvalues match...
1./eigA

# ... the eigenvalues of the inverse(A)
eigvals(inv(A))

## Symmetric positive definite matrix:
B = rand(3,3);
A = B'*B + 0.0001*eye(3)

# Has a matrix square root
sqrtA = sqrtm(A)

# Notice that the eigenvalues of sqrt(A)...
eigvals(sqrtA)

# Are equal to the sqrt of the eigs of A:
sqrt.(eigvals(A))

# Given A, find R s.t. A = R U:
A = rand(3,3)

# Find U first, then R:
U2 = A'*A;
U  = sqrtm(U2)  # sqrt root of U*U

# Now A=RU --> R=A invU
R = A*inv(U)  # not efficient algorithm (for now), but true

# Verify R is orthogonal:
R'*R

# Factorization A=RU is accurate:
norm(A-R*U)/norm(A)

?(svd)

A = rand(4,2)

# Thin SVD:
(U,S,V) = svd(A)

U

S

diagm(S)

V'

A

U*diagm(S)*V'

## Relative Error:
norm(U*diagm(S)*V'-A)/norm(A)



function vander(N)
    V = [(i/N)^(j-1) for i=1:N, j=1:N];
end

N = 10;
V = vander(N)

sv = svdvals(V)

## Matrix two-norm is largest singval:
norm(V)

sv[1]

cond(V)

## Definition of condition number in singular values:
sv[1]/sv[length(sv)]



using Plots
N = 1000;
V = vander(N) / norm(vander(N));  ## normalize V for clarity
sv = svdvals(V);
plot(log10.(sv))

A0  = V'*V;  
alpha = 1.e-10; 
A   = A0 + alpha*eye(N);
sv0 = svdvals(A0);
sv  = svdvals(A);
plot(log10.([sv0, sv]))

A = [V; sqrt(alpha)*eye(N)];

plot(log10.([svdvals(V), svdvals(A)]))

plot(log10.([sv0, svdvals(A'*A)]))