<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Linear Algebra Examples
This just shows the machanics of linear algebra calculations with python. See Lecture 5 for motivation and understanding.
Step1:  Resources
Tutorial for scipy.linalg
Exact solution of linear system of equations
\begin{align}
x + 2y &= 3 \
3x + 4y &= 17
\end{align}
Step2:  Using solve is faster and more stable numerically than using matrix inversion
Step3:  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
Step4:  Basic information about a matrix
Step5:  Trace
Step6:  Determinant
Step7:  Rank
Step8:  Norm
Step9:  Least-squares solution
Step10:  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
Step11:  Matrix Decompositions
Step12:  LU
Step13:  Choleskey
Step14:  QR
Step15:  Spectral
Step16:  Inverting A
Step17:  Powers of A
Step18:  SVD
<ASSISTANT_TASK:>
Python Code:
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
<END_TASK>