import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
%matplotlib inline

plt.style.use('ggplot')

"""
Explanation: Linear Algebra Examples
This just shows the machanics of linear algebra calculations with python. See Lecture 5 for motivation and understanding.
End of explanation
"""


A = np.array([[1,2],[3,4]])
A

b = np.array([3,17])
b

x = la.solve(A, b)
x

np.allclose(A @ x, b)

A1 = np.random.random((1000,1000))
b1 = np.random.random(1000)

"""
Explanation: Resources

Tutorial for scipy.linalg

Exact solution of linear system of equations
\begin{align}
x + 2y &= 3 \
3x + 4y &= 17
\end{align}
End of explanation
"""


%timeit la.solve(A1, b1)

%timeit la.inv(A1) @ b1

"""
Explanation: Using solve is faster and more stable numerically than using matrix inversion
End of explanation
"""


import scipy.linalg.lapack as lapack

lu, piv, x, info = lapack.dgesv(A, b)
x

"""
Explanation: Under the hood (Optional)
The solve function uses the dgesv fortran function to do the actual work. Here is an example of how to do this directly with the lapack function. There is rarely any reason to use blas or lapack functions directly becuase the linalg package provides more convenient functions that also perfrom error checking, but you can use Python to experiment with lapack or blas before using them in a language like C or Fortran.

How to interpret lapack function names
Summary of BLAS functions
Sumary of Lapack functions
End of explanation
"""


C = np.array([[1, 2+3j], [3-2j, 4]])
C

C.conjugate()

"""
Explanation: Basic information about a matrix
End of explanation
"""


def trace(M):
    return np.diag(M).sum()

trace(C)

np.allclose(trace(C), la.eigvals(C).sum())

"""
Explanation: Trace
End of explanation
"""


la.det(C)

"""
Explanation: Determinant
End of explanation
"""


np.linalg.matrix_rank(C)

"""
Explanation: Rank
End of explanation
"""


la.norm(C, None) # Frobenius (default)

la.norm(C, 2) # largest sinular value

la.norm(C, -2) # smallest singular value

la.svdvals(C)

"""
Explanation: Norm
End of explanation
"""


la.solve(A, b)

x, resid, rank, s = la.lstsq(A, b)
x

A1 = np.array([[1,2],[2,4]])
A1

b1 = np.array([3, 17])
b1

try:
    la.solve(A1, b1)
except la.LinAlgError as e:
    print(e)

x, resid, rank, s = la.lstsq(A1, b1)
x

A2 = np.random.random((10,3))
b2 = np.random.random(10)

try:
    la.solve(A2, b2)
except ValueError as e:
    print(e)

x, resid, rank, s = la.lstsq(A2, b2)
x

"""
Explanation: Least-squares solution
End of explanation
"""


def least_squares(X, y):
    return la.solve(X.T @ X, X.T @ y)

least_squares(A2, b2)

"""
Explanation: Normal equations
One way to solve least squares equations $X\beta = y$ for $\beta$ is by using the formula $\beta = (X^TX)^{-1}X^Ty$ as you may have learnt in statistical theory classes (or can derive yourself with a bit of calculus). This is implemented below.
Note: This is not how the la.lstsq function solves least square problems as it can be inefficent for large matrices.
End of explanation
"""


A = np.array([[1,0.6],[0.6,4]])
A

"""
Explanation: Matrix Decompositions
End of explanation
"""


p, l, u = la.lu(A)

p

l

u

np.allclose(p@l@u, A)

"""
Explanation: LU
End of explanation
"""


U = la.cholesky(A)
U

np.allclose(U.T @ U, A)

# If workiing wiht complex matrices
np.allclose(U.T.conj() @ U, A)

"""
Explanation: Choleskey
End of explanation
"""


Q, R = la.qr(A)

Q

np.allclose((la.norm(Q[:,0]), la.norm(Q[:,1])), (1,1))

np.allclose(Q@R, A)

"""
Explanation: QR
End of explanation
"""


u, v = la.eig(A)

u

v

np.allclose((la.norm(v[:,0]), la.norm(v[:,1])), (1,1))

np.allclose(v @ np.diag(u) @ v.T, A)

"""
Explanation: Spectral
End of explanation
"""


np.allclose(v @ np.diag(1/u) @ v.T, la.inv(A))

"""
Explanation: Inverting A
End of explanation
"""


np.allclose(v @ np.diag(u**5) @ v.T, np.linalg.matrix_power(A, 5))

"""
Explanation: Powers of A
End of explanation
"""


U, s, V = la.svd(A)

U

np.allclose((la.norm(U[:,0]), la.norm(U[:,1])), (1,1))

s

V

np.allclose((la.norm(V[:,0]), la.norm(V[:,1])), (1,1))

np.allclose(U @ np.diag(s) @ V, A)

"""
Explanation: SVD
End of explanation
"""


