import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline


M = np.mat("3 -2; 1 0")
w, V = np.linalg.eig(M)

w

V

"""
Explanation: Linear Algebra via Examples
We demonstrate a few important LA concepts via examples. 
Eigenvalues and Eigenvectors
End of explanation
"""


V[:, 0]

3 * V[:, 0]

for i in range(2):
    print('i = {}'.format(i+1))
    print('Left  = ', (M @ V[:, i]).reshape(1,2))
    print('Right = ', (w[i] * V[:, i]).reshape(1,2))

"""
Explanation: Verify the property that $$\mathbf{M} \mathbf{v} = \lambda \mathbf{v} $$
End of explanation
"""


U = np.mat("1 0; 0 1")
V = M @ U
V

U2 = np.hstack([np.array([i-1, j-1]).reshape(2,1) for i in range(3) for j in range(3)])
U2

V2 = M @ U2
V2

plt.plot(U2[0,:], U2[1,:], 'bx')
plt.title('original space')
plt.draw()
plt.show()

plt.plot(V2[0,:], V2[1,:], 'r+')
plt.title('new space')
plt.draw()
plt.show()

plt.plot(U2[0,:], U2[1,:], 'bx')
plt.plot(V2[0,:], V2[1,:], 'r+')
plt.title('plot together')
plt.draw()
plt.show()

"""
Explanation: Impact to basis vectors
End of explanation
"""


from scipy import linalg

A = np.array([ [5, 3], [0, -1], [12, 1], [4, 2] ])
A

"""
Explanation: With some more efforts, one can draw the original coordinate systems and the transformed one in a more visual way. This will hopefully strengthen the idea that a matrix multiplication represents a linear transformation.
SVD
End of explanation
"""


U, s, Vt = linalg.svd(A, full_matrices=False)
S = np.diag(s)
V = Vt.T

U

s

S

Vt

V

"""
Explanation: We only do a reduced SVD: $$\mathbf{A} = \mathbf{U} \mathbf{S} \mathbf{V}^{\top}, $$ 
so that $\mathbf{U}$ is of the same shape as $\mathbf{A}$.
End of explanation
"""


print(U.shape, S.shape, Vt.shape)

reconA = np.round(U @ S @ Vt, 5)
reconA

"""
Explanation: We demonstrate a few useful properties of SVD: 

(1) Exact Recovery
End of explanation
"""


U1 = U[:, 0]
U1

s1 = s[0]
s1

Vt1 = Vt[0, :]
Vt1

reconA1 = np.round(s1 * np.outer(U1, Vt1), 5)
reconA1

A

"""
Explanation: (2) Approximation

Another view of SVD is to decompose the original matrix into a sum of rank-1 matrices of decreasing "importance", i.e., 
$$
    \mathbf{A} = \sum_{i} {\sigma_i} \left( \mathbf{U}_i \otimes \mathbf{V}_i \right), 
$$
where $\otimes$ is the vector outer product, i.e., 
$$
\mathbf{u} \otimes \mathbf{v} = \mathbf{u} \otimes \mathbf{v}^{\top}
$$
In the following, we will approximate $\mathbf{A}$ by only one summand.
End of explanation
"""


a = np.random.rand(8)
b = np.random.rand(8)
sigma = 0.01
At = np.array( [a + sigma*np.random.randn(8) for _ in range(4)])
At = np.vstack((At, b))
A = At.T
print(a)
print(b)
print()
print(A)


def approximate_one(A):
    U, s, Vt = linalg.svd(A, full_matrices=False)
    print('singular values are: \n', s) 
    S = np.diag(s)
    V = Vt.T
    # only 1 summand
    U1 = U[:, 0]
    s1 = s[0]
    Vt1 = Vt[0, :]
    A1 = s1 * np.outer(U1, Vt1)
    return A1

A1 = approximate_one(A)
A1

print(A1 - A)

[ np.sum( (A1-A)[:, i] ) for i in range(A.shape[1]) ]

"""
Explanation: (3) Feature reduction

We manually creat a matrix $\mathbf{A}$ such that all but its last column is heavily correlated.
End of explanation
"""


