#!/usr/bin/env python
# coding: utf-8

# # 线性代数模块
# 
# scipy.linalg是SciPy中负责线性代数计算的模块。NumPy也有numpy.linalg模块，不过建议使用scipy.linalg，原因有二：
# - scipy.linalg包含numpy.linalg中的所有函数，同时还包含了很多numpy.linalg中没有的函数；
# - scipy.linalg默认使用线性代数库BLAS/LAPACK等对运算进行加速，而在numpy.linalg中，这些加速是需要自己配置，不是默认设置:

# In[1]:


import numpy.linalg


# In[2]:


import scipy.linalg


# In[3]:


len(dir(numpy.linalg))


# In[4]:


len(dir(scipy.linalg))


# In[5]:


import numpy as np


# In[6]:


from scipy import linalg


# ## 基本的矩阵操作

# 在NumPy中，矩阵有两种表示方法：矩阵类型和二维数组类型。
# 
# 矩阵类型可以用`np.mat()`或者`np.matrix()`创建：

# In[7]:


A = np.mat("[1, 2; 3, 4]")


# In[8]:


A


# In[9]:


A = np.matrix("[1, 2; 3, 4]")


# In[10]:


A


# 转置矩阵为：

# In[11]:


A.T


# 矩阵的逆矩阵：

# In[12]:


A.I


# 矩阵乘法：

# In[13]:


A.I * A


# In[14]:


A * A.I


# 数组类型也可以完成类似的操作：

# In[15]:


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


# In[16]:


A


# 转置：

# In[17]:


A.T


# 其逆矩阵可以用`linalg.inv()`函数计算：

# In[18]:


A.dot(linalg.inv(A))


# In[19]:


A @ linalg.inv(A)


# In[20]:


linalg.inv(A) @ A


# `scipy.linalg`也可以作用在矩阵类型上，返回的还是一个数组：

# In[21]:


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


# In[22]:


linalg.inv(A)


# 矩阵类型的使用比较少，通常都用数组类型替代矩阵类型进行计算。

# ## 线性方程组的求解

# 考虑下列方程组：
# ```
#  x + 3y + 5z = 10
# 2x + 5y -  z = 8
# 2x + 3y + 8z = 3
# ```
# 
# 矩阵形式为：
# $$
# \left[
# \begin{array}{1}
# 1 & 3 & 5 \\
# 2 & 5 & -1 \\
# 2 & 3 & 8 \\
# \end{array}
# \right]
# \left[
# \begin{array}{1}
# x \\
# y \\
# z \\
# \end{array}
# \right]
# =
# \left[
# \begin{array}{1}
# 10 \\
# 8 \\
# 3 \\
# \end{array}
# \right]
# $$
# 
# 其解为：
# $$
# \left[
# \begin{array}{1}
# x \\
# y \\
# z \\
# \end{array}
# \right]
# =
# \left[
# \begin{array}{1}
# 1 & 3 & 5 \\
# 2 & 5 & -1 \\
# 2 & 3 & 8 \\
# \end{array}
# \right]^{-1}
# \left[
# \begin{array}{1}
# 10 \\
# 8 \\
# 3 \\
# \end{array}
# \right]
# $$
# 
# 故：

# In[23]:


A = np.array([[1, 3, 5],
              [2, 5, -1],
              [2, 3, 8]])


# In[24]:


b = np.array([10, 8, 3])


# 解：

# In[25]:


linalg.inv(A) @ b


# 也可以直接使用`linalg.solve()`函数，直接求解这个方程：

# In[26]:


linalg.solve(A, b)


# ## 行列式的计算：

# `linalg.det()`函数可以计算矩阵的行列式：

# In[27]:


linalg.det(A)


# 行列式非零表示矩阵是可逆的。逆矩阵的行列式与矩阵的行列式满足相乘等于1的关系：

# In[28]:


linalg.det(linalg.inv(A))


# In[29]:


linalg.det(linalg.inv(A)) * linalg.det(A)


# ## 范数的计算

# In[30]:


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


# 默认情况下，`linalg.norm()`函数计算的是矩阵的Frobenius范数，该范数计算的是矩阵所有元素平方和的平方根：

# In[31]:


linalg.norm(A)


# In[32]:


linalg.norm(A, "fro")


# 可以通过`linalg.norm()`函数中的第二个参数，来修改矩阵范数的计算方法。例如，矩阵的1范数，计算的是矩阵中每列元素的和的最大值：

# In[33]:


linalg.norm(A, 1)


# 矩阵的-1范数，计算的是矩阵中每列元素的和的最小值：

# In[34]:


linalg.norm(A, -1)


# 矩阵的2范数，计算的是矩阵的最大奇异值：

# In[35]:


linalg.norm(A, 2)


# 矩阵的-2范数，计算的是矩阵的最小奇异值：

# In[36]:


linalg.norm(A, -2)


# 矩阵的正无穷范数，计算的是矩阵中每行元素的和的最大值：

# In[37]:


linalg.norm(A, np.inf)


# 矩阵的负无穷范数，计算的是矩阵中每行元素的和的最小值：

# In[38]:


linalg.norm(A, -np.inf)


# `linalg.norm()`函数也可以用来计算向量的范数（也称模）：

# In[39]:


b = np.array([10, 8, 3])


# 对于向量来说，linalg.norm()函数默认是向量的2范数，即计算向量元素的平方和的平方根：

# In[40]:


linalg.norm(b)


# In[41]:


linalg.norm(b, 2)


# ## 特征值分解

# 对于方阵A，特征值λ和特征向量v满足：`Av=λv`。函数`linalg.eig()`可以用来求解特征值和特征向量：

# In[42]:


A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])


# In[43]:


l, v = linalg.eig(A)


# In[44]:


np.allclose(A.dot(v), l * v)


# ## 广义逆
# 
# 对于形状为m×n的矩阵A，可以用`linalg.pinv()`求解其广义逆。A的广义逆B满足：`A=ABA`，`B=ABA`：

# In[45]:


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


# In[46]:


B = linalg.pinv(A)


# In[47]:


np.allclose(A, A @ B @ A)


# In[48]:


np.allclose(B, B @ A @ B)


# ## 奇异值分解

# 对于形状为m×n的矩阵A，其奇异值分解满足$A=U\SigmaV^H$，矩阵Σ形状为M×N，主对角线上的元素被称为奇异值。U、V分别为M阶、N阶正交单位矩阵。
# 
# 函数`linalg.svd()`可以对矩阵进行奇异值分解：

# In[49]:


U, s, Vh = linalg.svd(A)


# In[50]:


U


# In[51]:


s


# In[52]:


Vh


# 从奇异值恢复矩阵$\Sigma$：

# In[53]:


S = linalg.diagsvd(s, 2, 3)


# In[54]:


np.allclose(A, U @ S @ Vh)


# 正交矩阵：

# In[55]:


U @ U.T


# In[56]:


Vh @ Vh.T


# In[ ]: