%% Eigenvalues and Singular Values
% Symmetric matrices have real eigenvalues and a set of orthonormal eigenvectors. What can be said about matrics that are not symmetric?
% First look at the EVD.

clear
A=gallery(3)
[V,D] = eig(A)
norm(A-V*D*V')

%
% Now look at the SVD.

[U,S,V] = svd(A)
norm(A - U*S*V')

%
% Let's create a symmetric matrix from a known set of eigenvalues.

n = 3;
lambda = (1:n)';
D = diag(lambda);
[V,~] = qr(randn(n));
A = V*D*V'

%
% We can perform an SVD and will see that U=V.

[U,S,V] = svd(A)

%
% The Rayleigh quotient can be used to get good estimates of eigenvalues. Recall that R_A(\mathbf{x})=\frac{\mathbf{x}^*A\mathbf{x}}{\mathbf{x}^*\mathbf{x}}. Here we apply it to an exact eigenvector of A.

R = @(x) (x'*A*x)/(x'*x);
R(V(:,2))

%
% Now, perturb the eigenvector and see how the eigenvalue is affected.

epsilon=1e-2;
R(V(:,2)+randn(n,1)*epsilon)

%
% Now we consider the geometric transfromations of dialation (stretching/shrinking) and rotation on unit vectors. First we show dilations.

clear                      % Get vectors |x|=1.
t=linspace(0,2*pi,100);
x1=cos(t);
x2=sin(t);
x=[x1;x2];

A1=[1 0;0 2];              % Stretching
ys = A1*x;
figure
subplot(1,2,1) 
plot(x1,x2); axis equal
subplot(1,2,2) 
plot(ys(1,:), ys(2,:)); axis equal

%
% Next, we rotate the vectors using the 2D rotation matrix. R=\left( \begin{array}{cc} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta\end{array}  \right).

th=pi/6;                                   % Rotation angle
R=[cos(th) -sin(th); sin(th) cos(th)];
A=R*A1
[V,D]=eig(A)
y = A*x;
figure
p(1)=subplot(1,2,1) 
plot(ys(1,:), ys(2,:)); %axis equal
p(2)=subplot(1,2,2) 
plot(y(1,:), y(2,:)); %axis equal
linkaxes(p,'y')
ax = axis;
axis([ax(1:2) min(ys(2,:)), max(ys(2,:))])

%
% Plot eigenvectors for a symmetric matrix.

A=[3 1;1 2];   % Symmetric
[V,D]=eig(A)
y = A*x;
figure
plot(y(1,:), y(2,:)); axis equal
hold on
line ([0,-V(1,1)],[0,-V(2,1)]);
line ([0,-V(1,2)],[0,-V(2,2)]);
hold off

%
% Note that the eigenvectors are orthonormal.

V(:,1)'*V(:,2)
V(:,1)'*V(:,1)
V(:,2)'*V(:,2)

%
% What happens in the non-symmetric case? The eigenvectors do not point along directions of stretching and shrinking.

A=[1 2;3 4];   % Not Symmetric
A=[3 2;0 2];
%A=[3 1;6 2];
y = A*x;
[V,D]=eig(A)

figure
plot(y(1,:), y(2,:)); axis equal
hold on
line ([0,V(1,1)],[0,V(2,1)]);
line ([0,V(1,2)],[0,V(2,2)]);
hold off

V^(-1)*[2;2]

%
% Here we note that A^*Ais symmetric. So, it has real eigenvalues and orthogonasl eigenvectors. Are they useful?

y2 = A'*A*x;
[V2,D2]=eig(A'*A)

%figure
plot(y2(1,:), y2(2,:)); axis equal
hold on
plot([0,-D(2,2)*V2(1,1)],[0,-D(2,2)*V2(2,1)],'r');
plot([0,-D(1,1)*V2(1,2)],[0,-D(1,1)*V2(2,2)],'r');
hold off

%

figure
plot(y(1,:), y(2,:)); axis equal
hold on
line ([0,V(1,1)],[0,V(2,1)]);
line ([0,V(1,2)],[0,V(2,2)]);
plot([0,-D(2,2)*V2(1,1)],[0,-D(2,2)*V2(2,1)],'r');
plot([0,-D(1,1)*V2(1,2)],[0,-D(1,1)*V2(2,2)],'r');
hold off
V2(1,1)^2+V2(2,1)^2

%
% One can do a similar study in 3D. The set of unit vectors lie on a unit sphere and dilations and rotations transform the sphere into an ellipsoid.

figure
[X,Y,Z] = sphere;
surf(X,Y,Z)
axis equal

%
% We create a symmetric matrix.

n = 3;
lambda = (1:n)';
D = diag(lambda);
[V,~] = qr(randn(n));
A = V*D*V'
[V,D]=eig(A)

%
% It is interesting to use the quadratic form in three variables. We want to plot ax^2+by^2+cz^2+dxy+exz+fyz+gx+hy+jz=r^2. So, we not the form
%   \left( \begin{array}{cccc} x& y& z& 1 \end{array}\right) \left( \begin{array}{cccc} a  & \frac d2 & \frac e2 & \frac g2 \\ \frac d2  & b & \frac f2 & \frac h2 \\ \frac e2  & \frac f2 & c & \frac j2 \\ \frac g2  & \frac h2 & \frac j2 & 1\ \end{array} \right)  \left( \begin{array}{c} x \\ y\\
%  z\\ 1 \end{array}\right) = r^2+1.
% Let  \mathbf{x}= \left( \begin{array}{cccc} x,& y,& z,& 1 \end{array}\right) ^T and M=\left( \begin{array}{cccc} a  & d & e & g \\ 0  & b & f & h \\ 0  & 0 & c & j \\ 0  & 0& 0 & 1\ \end{array} \right).Then, this can be written more compactly as \frac 12\mathbf{x}^T(M+M^T)\mathbf{x}=r^2+1.
% The code for doing this is in the function below. Let's say we want to plot the ellipsoid x^2+2y^2+3z^2+4xy+5xz+6yz+7x+8y+9z=10.The ellipsoid function is called as follows:

a=1;b=2;c=3;d=4;e=5;f=6;
g=7;h=8;j=9;R=sqrt(10);
ellipsoid(a, b, c, d, e, f, g, h, j, R)

%
% Next we consider no translation and write Q=\frac 12(M+M^T)  from the EVD as the average of the above real upper triangular matrix and its transpose. Then , we can introduce the appropriate parameters to study \mathbf{x}^TQ\mathbf{x}=r^2+1.
% [ Need to relate Q to the dilation+rotation matrix vs V. ]

a=V(1,1);  b=V(2,2);  c=V(3,3);
d=V(1,2);  e=V(1,3);  f=V(2,3);
g=0;  h=0;  j=0;  R=sqrt(1);
ellipsoid(a, b, c, d, e, f, g, h, j, R)

% Angle of rotation: 1+2 cos(theta) = tr(V)
acos((trace(V)-1)/2)*180/pi
det(V)             % Not a proper rotation
A*V(:,1)-V(:,1)    % Axis of rotation
hold on
k=3;
line ([-k*V(1,1),k*V(1,1)],[-k*V(2,1),k*V(2,1)],[-k*V(3,1),k*V(3,1)],'LineWidth',3,'Color','b');
hold off

%
% The end so far.

function ellipsoid(a, b, c, d, e, f, g, h, j, R)
M = (1/R^2) * [a d e g;0 b f h; 0 0 c j; 0 0 0 1];
M = (M + M')/2;
S = abs(M).^0.5;

n = 30;
[th, ph] = meshgrid(linspace(0, 2*pi, n), linspace(0, pi, n));
spherepoints = zeros(4, n, n);
spherepoints(1, :,:) = cos(th) .* sin(ph);
spherepoints(2, :,:) = sin(th) .* sin(ph);
spherepoints(3, :,:) = cos(ph);
spherepoints(4, :, :) = 1;
sp = reshape(spherepoints, 4, []); 
epoints = S \ sp;
epoints = reshape(epoints, 4, n, n);

surf(squeeze(epoints(1, :,:)),...
    squeeze(epoints(2,:,:)),...
    squeeze(epoints(3,:,:)));
set(gca, 'DataAspectRatio', [1,1,1]);
figure(gcf);
end