Linear Algebra Examples
====

This just shows the machanics of linear algebra calculations with python. See Lecture 5 for motivation and understanding.


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


```python
plt.style.use('ggplot')
```

Resources
----

- [Tutorial for `scipy.linalg`](http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html)

Exact solution of linear system of equations
----

\begin{align}
x + 2y &= 3 \\
3x + 4y &= 17
\end{align}



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




    array([[1, 2],
           [3, 4]])




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




    array([ 3, 17])




```python
x = la.solve(A, b)
x
```




    array([11., -4.])




```python
np.allclose(A @ x, b)
```




    True




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

### Using solve is faster and more stable numerically than using matrix inversion


```python
%timeit la.solve(A1, b1)
```

    437 ms ± 115 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)



```python
%timeit la.inv(A1) @ b1
```

    584 ms ± 127 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


### 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](http://www.netlib.org/lapack/lug/node24.html)
- [Summary of BLAS functions](http://cvxopt.org/userguide/blas.html)
- [Sumary of Lapack functions](http://cvxopt.org/userguide/lapack.html)


```python
import scipy.linalg.lapack as lapack
```


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

Basic information about a matrix
----


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




    array([[1.+0.j, 2.+3.j],
           [3.-2.j, 4.+0.j]])




```python
C.conjugate()
```




    array([[1.-0.j, 2.-3.j],
           [3.+2.j, 4.-0.j]])



#### Trace


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


```python
trace(C)
```




    (5+0j)




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




    True



#### Determinant


```python
la.det(C)
```




    (-8-5j)



#### Rank


```python
np.linalg.matrix_rank(C)
```




    2



#### Norm


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




    6.557438524302




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




    6.389028023601217




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




    1.4765909770949925




```python
la.svdvals(C)
```

Least-squares solution
----


```python
la.solve(A, b)
```


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


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


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


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


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


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


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


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

### 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.


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


```python
least_squares(A2, b2)
```

Matrix Decompositions
----


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




    array([[1. , 0.6],
           [0.6, 4. ]])



### LU


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


```python
p
```




    array([[1., 0.],
           [0., 1.]])




```python
l
```




    array([[1. , 0. ],
           [0.6, 1. ]])




```python
u
```




    array([[1.  , 0.6 ],
           [0.  , 3.64]])




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




    True



### Choleskey


```python
U = la.cholesky(A)
U
```




    array([[1.       , 0.6      ],
           [0.       , 1.9078784]])




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




    True




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




    True



### QR


```python
Q, R = la.qr(A)
```


```python
Q
```




    array([[-0.85749293, -0.51449576],
           [-0.51449576,  0.85749293]])




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




    True




```python

```


```python
np.allclose(Q@R, A)
```




    True



### Spectral

When matrix is symmetric, you can use la.eigh


```python
u, v = la.eig(A)
```


```python
u
```




    array([0.88445056+0.j, 4.11554944+0.j])




```python
v
```




    array([[-0.98195639, -0.18910752],
           [ 0.18910752, -0.98195639]])




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




    True




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




    True



#### Inverting A


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

#### Powers of A


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

### SVD


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


```python
U
```


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


```python
s
```


```python
V
```


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


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


```python

```

2a1d4d98a850642594e1987d5c8c14440dae0af8
ipynb
Jupyter Notebook
notebooks/copies/lectures/T06_Linear_Algebra_Examples.ipynb
robkravec/sta-663-2021
4dc8018f7b172eaf81da9edc33174768ff157939
MIT
notebooks/copies/lectures/T06_Linear_Algebra_Examples.ipynb
robkravec/sta-663-2021
4dc8018f7b172eaf81da9edc33174768ff157939
MIT
notebooks/copies/lectures/T06_Linear_Algebra_Examples.ipynb
robkravec/sta-663-2021
4dc8018f7b172eaf81da9edc33174768ff157939
MIT
Qwen/Qwen-72B
1. YES
2. YES