Skip to content

Latest commit

 

History

History
921 lines (475 loc) · 19.3 KB

线性代数.md

File metadata and controls

921 lines (475 loc) · 19.3 KB
Raw

线性代数

numpy中的多维数组最常见的用法就是解决线性代数问题,尤其是矩阵运算.

向量间计算

向量间计算多数都是对应项操作,这个ufunc可以直接支持.比较特殊的操作是点乘和求夹角

import numpy as np

求模

向量中模的定义为各个项的平方和开方

$$ |v| = \sqrt {\sum_{i=0}^{N} {v_i^2}} $$

numpy中使用np.linalg.norm(v)来计算

v_1 = np.array([3,5])
l_v1 = np.linalg.norm(v_1)
print(l_v1)
5.830951894845301

点积(内积)

向量的点积为各项相乘的和,为一个标量.numpy中使用np.linalg.dot来计算其值,或者使用python3.6定义的矩阵乘法符号@来计算

$$ a \cdot b = \sum_{i=0}^{n} a_i\cdot b_i $$

v_1 = np.array([3,5])
v_2 = np.array([4,2])
v_1@v_2
22

叉乘(外积)

$$ |C| = | v_1 \times v_2 |= |v_1| |v_2|sin<a,b> $$

或者在二维情况下可以看做:

$$ (x_1,y_1) \times (x_2,y_2)= x_1y_2 -x_2y_1 $$

C的方向用右手定则,右手4指从v_1不超过180度转向v_2时大拇指的方向即为C的方向.其值为新向量的模

numpy中没有直接计算的方法,这边需要迂回下,通过计算余弦,再转换为角度计算正弦再求出

v_1 = np.array([3,5])
v_2 = np.array([4,2])

v_1_norm = np.linalg.norm(v_1)
v_2_norm = np.linalg.norm(v_2)
cos_theta = (v_1@v_2)/(v_1_norm * v_2_norm) # 计算余弦
theta = np.arccos(cos_theta)
cross_product = np.sin(theta)* (v_1_norm * v_2_norm)
cross_product
14.000000000000002

矩阵计算

通常我们认为向量是一维的矩阵,numpy中通常使用二维数组代表矩阵,概括下矩阵计算的工具包括

向量,矩阵运算:

函数 说明
dot 向量乘法
vdot 向量点乘
inner 向量内积
outer 向量外积
matmul 矩阵乘法
trnsordot 张量乘法
einsum 评估操作数上的爱因斯坦求和约定
linalg.matrix_power 矩阵幂
kron 克罗内克积

矩阵分解

函数 说明
linalg.cholesky(a) Cholesky 分解
linalg.qr(a[, mode]) QR分解
linalg.svd(a[, full_matrices, compute_uv]) 奇异值分解

矩阵征值操作

函数 说明
linalg.norm(x[, ord, axis, keepdims]) 矩阵或向量范数
linalg.cond(x[, p]) 计算矩阵的条件数
linalg.det(a) 计算矩阵行列式
linalg.matrix_rank(M[, tol]) 使用SVD方法返回阵列的矩阵秩
linalg.slogdet(a) 计算数组行列式的符号和（自然）对数
trace(a[, offset, axis1, axis2, dtype, out]) 计算对角线元素的和
diag 以一维数组的形式返回方阵的对角线(或非对角线)元素,或将一维数组转换为方阵(非对角线元素为0)
eig 计算方阵的本征值和本征向量

求解方程和求逆矩阵

函数 说明
linalg.solve(a, b) 解线性方程组Ax=b
linalg.tensorsolve(a, b[, axes]) 解张量表达式Ax = b
linalg.lstsq(a, b[, rcond]) 计算Ax=b的最小二乘解
linalg.inv(a) 计算方阵的逆
linalg.pinv(a[, rcond]) 计算矩阵的Moore-Penrose伪逆
linalg.tensorinv(a[, ind]) 计算N维数组的“逆”。 
A = np.arange(9).reshape(3,3)
print(A)
[[0 1 2]
 [3 4 5]
 [6 7 8]]

B = np.arange(13,22).reshape(3,3)
print(B)
[[13 14 15]
 [16 17 18]
 [19 20 21]]

矩阵的秩

线性代数中 $$ \mathbf {A} ={\begin{bmatrix}1&2\3&4\end{bmatrix}} $$

一个矩阵A的列秩是A的线性独立的纵列的极大数目.类似地,行秩是A的线性独立的横行的极大数目.

矩阵的列秩和行秩总是相等的,因此它们可以简单地称作矩阵A的秩.通常表示为r(A)，rk(A)或rank A.

矩阵的行秩与列秩相等,是线性代数基本定理的重要组成部分.其基本证明思路是,矩阵可以看作线性映射的变换矩阵,列秩为像空间的维度,行秩为非零原像空间的维度,因此列秩与行秩相等,即像空间的维度与非零原像空间的维度相等(这里的非零原像空间是指约去了零空间后的商空间:原像空间).这从矩阵的奇异值分解就可以看出来.

给出这一结果的两种证明.第一个证明是简短的,仅用到向量的线性组合的基本性质.第二个证明利用了正交性.第一个证明利用了列空间的基,第二个证明利用了行向量空间的基.第一个证明适用于定义在标量域上的矩阵,第二个证明适用于内积空间.二者都适用于实或复的欧氏空间,也都易于修改去证明当A是线性变换的情形.

M_1 = np.matrix(np.arange(1,5).reshape(2,2))
M_1
matrix([[1, 2],
        [3, 4]])

np.linalg.matrix_rank(M_1)
2

方阵的迹

迹就是方阵主对角线元素之和

np.trace(M_1)
5

转置矩阵(transpose)

将矩阵延对角线翻转

M_1.T
matrix([[1, 3],
        [2, 4]])

共轭矩阵(hermitian)(复数为元素)

M_2 = np.matrix([[1+1j,2-4j],[3-1j,2+3j]])
M_2
matrix([[1.+1.j, 2.-4.j],
        [3.-1.j, 2.+3.j]])

M_2.H
matrix([[1.-1.j, 3.+1.j],
        [2.+4.j, 2.-3.j]])

逆矩阵(inverse)

在线性代数中,给定一个n阶方阵$\mathbf{A}$，若存在一n阶方阵$ \mathbf {B}$,使得 $ \mathbf{AB}=\mathbf{BA}=\mathbf{I}_n$,其中 $ \mathbf{I}_n $为n阶单位矩阵,则称 $\mathbf{A} $是可逆的,且 $\mathbf {B} $是$\mathbf{A}$的逆矩阵,记作$\mathbf {A} ^{-1}$

只有正方形（n×n）的矩阵,即方阵,才可能但非必然有逆矩阵.若方阵 $\mathbf{A}$的逆矩阵存在,则称 $\mathbf{A}$为非奇异方阵或可逆方阵.

性质有:

  • $ \left (A^{-1} \right )^{-1}=A $
  • $ (\lambda A)^{-1}=\frac{1}{\lambda}\times A^{-1} $
  • $ (AB)^{-1}=B^{-1}A^{-1} $
  • $ \left (A^\mathrm{T} \right )^{-1}=\left (A^{-1} \right )^{\mathrm{T}}$ ($ A^{\mathrm{T}}$ 为A的转置)
  • $ \det(A^{-1})=\frac{1}{\det(A)} $（det为行列式）
M_1.I
matrix([[-2. ,  1. ],
        [ 1.5, -0.5]])

伴随矩阵(adjoint)

在线性代数中,一个方形矩阵的伴随矩阵是一个类似于逆矩阵的概念.如果矩阵可逆,那么它的逆矩阵和它的伴随矩阵之间只差一个系数.然而伴随矩阵对不可逆的矩阵也有定义,并且不需要用到除法.

对n×n的矩阵A和B,有:

  • $\mathrm{adj}(\mathbf{I}) = \mathbf{I}$

  • $(\mathbf{AB}) = \mathrm{adj}(\mathbf{B}),\mathrm{adj}(\mathbf{A})$

  • $\mathrm{adj}(\mathbf{A}^T) = \mathrm{adj}(\mathbf{A})^T $

  • $ \det\big(\mathrm{adj}(\mathbf{A})\big) = \det(\mathbf{A})^{n-1}$

  • $ \mathrm{adj}(k \mathbf{A}) = k^{n-1} \ \mathrm{adj}(\mathbf{A}) $

  • 当n>2时,$\mathrm{adj}(\mathrm{adj}(\mathbf{A})) =(\det \mathbf{A})^{n-2} \mathbf{A} $

  • 如果A可逆,那么 $ \mathrm{adj}(\mathbf{A}^{-1}) = \mathrm{adj}(\mathbf{A})^{-1} = \frac{A}{\det A} $

  • 如果A是对称矩阵,那么其伴随矩阵也是对称矩阵;如果A是反对称矩阵,那么当n为偶数时,A的伴随矩阵也是反对称矩阵,n为奇数时则是对称矩阵.

  • 如果A是(半)正定矩阵,那么其伴随矩阵也是(半)正定矩阵.

  • 如果矩阵A和B相似,那么$\mathrm{adj}(\mathbf{A})$和 $\mathrm{adj}(\mathbf{B})$也相似.

  • 如果n>2,那么非零矩阵A是正交矩阵当且仅当 $ \mathrm{adj}(\mathbf{A}) = \pm A^T $

  • 伴随矩阵的秩

    当矩阵A可逆时,它的伴随矩阵也可逆,因此两者的秩一样,都是n.当矩阵A不可逆时,A的伴随矩阵的秩通常并不与A相同.当A的秩为n-1时,其伴随矩阵的秩为1,当A的秩小于n-1时,其伴随矩阵为零矩阵.

  • 伴随矩阵的特征值

    设矩阵A在复域中的特征值为 $\lambda_1, \lambda_2 \cdots \lambda_n$ (即为特征多项式的n个根),则A的伴随矩阵的特征值为$\lambda_2 \lambda_3 \cdots \lambda_{n}, \ \lambda_1 \lambda_3 \cdots \lambda_{n}, \cdots , \lambda_1 \lambda_2 \cdots \lambda_{n-1} $

np.dot(np.linalg.det(M_1),M_1.I)
matrix([[ 4., -2.],
        [-3.,  1.]])

矩阵的范数(matrix norms)

范数(norm),是具有"长度"概念的函数.在线性代数中其为向量空间内的所有向量赋予非零的正长度或大小.半范数反而可以为非零的向量赋予零长度.

举一个简单的例子,一个二维度的欧氏几何空间 $\mathbb {R} ^{2}$就有欧氏范数.在这个向量空间的元素(譬如：(3,7))常常在笛卡儿坐标系统被画成一个从原点出发的箭号.每一个向量的欧氏范数就是箭号的长度.

拥有范数的向量空间就是赋范向量空间.同样,拥有半范数的向量空间就是赋半范向量空间.

np.linalg.norm(M_1)
5.477225575051661

和和差

矩阵计算和和差要求两个矩阵形状一致,就是对应下标的各项计算的结果,numpy的ufunc刚好可以满足

A+B
array([[13, 15, 17],
       [19, 21, 23],
       [25, 27, 29]])

A-B
array([[-13, -13, -13],
       [-13, -13, -13],
       [-13, -13, -13]])

矩阵与标量的积

矩阵与标量的积就是标量在各项上的积.这也刚好ufunc就可以直接支持.

A*3
array([[ 0,  3,  6],
       [ 9, 12, 15],
       [18, 21, 24]])

3*A
array([[ 0,  3,  6],
       [ 9, 12, 15],
       [18, 21, 24]])

矩阵的内积(点积)

矩阵的乘法必须前一个矩阵的行数与后一个举证的列数相同.返回的是一个以前一个矩阵的行数为行数,后一个矩阵的列数为列数的新矩阵.

即形状上 $ (x,n)@ (n,y) -> (x,y) $,n是一致的

C = np.arange(12).reshape(4,3)
C
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11]])

D = np.arange(1,7).reshape(3,2)
D
array([[1, 2],
       [3, 4],
       [5, 6]])

C@D
array([[ 13,  16],
       [ 40,  52],
       [ 67,  88],
       [ 94, 124]])

求特征值特征向量

M_lamida=np.matrix([[3,0,-1],[2,4,2],[-1,0,3]])

np.linalg.eig(M_lamida)
(array([4., 4., 2.]), matrix([[ 0.        ,  0.70710678,  0.40824829],
         [ 1.        ,  0.        , -0.81649658],
         [ 0.        , -0.70710678,  0.40824829]]))

第一项是特征值,第二项是特征向量

判断正定矩阵

正定矩阵的定义是:设M是n阶方阵,如果对任何非零向量z,都有$ z^TMz > 0$,其中$z^T$ 表示z的转置,就称M正定矩阵.

M_4=np.arange(16).reshape(4,4)
M_4
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])

M_4 = M_4+M_4.T#将方阵转换成对称阵
M_4
array([[ 0,  5, 10, 15],
       [ 5, 10, 15, 20],
       [10, 15, 20, 25],
       [15, 20, 25, 30]])

lambdas,_ = np.linalg.eig(M_4)
lambdas
array([ 6.74165739e+01, -7.41657387e+00,  1.82694656e-15, -1.72637110e-15])

#判断是否所有特征值都大于0
True if np.all(lambdas > 0) else False
False

因此矩阵不是正定矩阵

还有一种方式是使用cholesky分解的方法:

Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解.它要求矩阵的所有特征值必须大于零,故分解的下三角的对角元也是大于零的.

np.linalg.cholesky(np.arange(16).reshape(4,4))
---------------------------------------------------------------------------

LinAlgError                               Traceback (most recent call last)

<ipython-input-30-9d406430339a> in <module>
----> 1 np.linalg.cholesky(np.arange(16).reshape(4,4))


~/Lib/conda/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py in cholesky(a)
    757     t, result_t = _commonType(a)
    758     signature = 'D->D' if isComplexType(t) else 'd->d'
--> 759     r = gufunc(a, signature=signature, extobj=extobj)
    760     return wrap(r.astype(result_t, copy=False))
    761 


~/Lib/conda/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_nonposdef(err, flag)
     98 
     99 def _raise_linalgerror_nonposdef(err, flag):
--> 100     raise LinAlgError("Matrix is not positive definite")
    101 
    102 def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):


LinAlgError: Matrix is not positive definite

np.linalg.cholesky(M_4)
---------------------------------------------------------------------------

LinAlgError                               Traceback (most recent call last)

<ipython-input-31-d7b2b7d13630> in <module>
----> 1 np.linalg.cholesky(M_4)


~/Lib/conda/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py in cholesky(a)
    757     t, result_t = _commonType(a)
    758     signature = 'D->D' if isComplexType(t) else 'd->d'
--> 759     r = gufunc(a, signature=signature, extobj=extobj)
    760     return wrap(r.astype(result_t, copy=False))
    761 


~/Lib/conda/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_nonposdef(err, flag)
     98 
     99 def _raise_linalgerror_nonposdef(err, flag):
--> 100     raise LinAlgError("Matrix is not positive definite")
    101 
    102 def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):


LinAlgError: Matrix is not positive definite

报错了,因此可以看出不是正定的

我们试试测试一个单位矩阵

i=np.eye(4)
np.linalg.cholesky(i)
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

奇异值分解

arr = np.arange(9).reshape((3, 3)) + np.diag([1, 0, 1])
arr
array([[1, 1, 2],
       [3, 4, 5],
       [6, 7, 9]])

uarr, spec, vharr = np.linalg.svd(arr)
uarr
array([[-0.1617463 , -0.98659196,  0.02178164],
       [-0.47456365,  0.09711667,  0.87484724],
       [-0.86523261,  0.13116653, -0.48390895]])

spec
array([14.88982544,  0.45294236,  0.29654967])

vharr
array([[-0.45513179, -0.54511245, -0.70406496],
       [ 0.20258033,  0.70658087, -0.67801525],
       [-0.86707339,  0.45121601,  0.21115836]])

矩阵QR分解

X = np.random.randn(5,5)
X
array([[-0.18551814, -0.20917874, -0.40211628, -1.19894807, -1.46630053],
       [ 1.29711693,  0.20234821,  0.32854939, -1.8502084 , -1.30213228],
       [-0.89666881, -1.69630947, -0.0259825 , -1.45345265,  0.06733376],
       [-0.43573678, -1.83826832,  1.7405289 ,  0.49307392,  1.21645164],
       [-0.34511311, -0.48444181, -1.52753147,  0.20559079, -1.36148696]])

mat = X.T.dot(X)
mat
array([[ 2.82991386,  2.79049186,  0.29282322, -1.16004689, -1.53755427],
       [ 2.79049186,  6.57608063, -2.26488911,  1.33591464, -1.64758744],
       [ 0.29282322, -2.26488911,  5.63311054,  0.45615902,  4.35704246],
       [-1.16004689,  1.33591464,  0.45615902,  7.25866167,  4.38925923],
       [-1.53755427, -1.64758744,  4.35704246,  4.38925923,  7.18352088]])

q,r=np.linalg.qr(mat)
q
array([[-0.63936271,  0.21100986, -0.18461235, -0.49914373,  0.51328742],
       [-0.63045609, -0.56252424, -0.23227077,  0.4278655 , -0.22152427],
       [-0.06615758,  0.53636822, -0.63382613, -0.02263921, -0.55288716],
       [ 0.2620895 , -0.58615552, -0.32811236, -0.66437978, -0.19665346],
       [ 0.34737978, -0.08815989, -0.63448821,  0.35478672,  0.5857527 ]])

r
array([[-4.42614784, -6.0024363 ,  2.50112288,  3.29643578,  5.37932372],
       [ 0.        , -4.96300058,  3.70576854, -5.3932577 , -0.26673908],
       [ 0.        ,  0.        , -6.01256737, -5.55185025, -8.09309891],
       [ 0.        ,  0.        ,  0.        , -2.12496231, -0.40364251],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.51142418]])

张量计算

将数据继续扩展,那么矩阵,向量,常数就都是张量的特殊形式了.

张量(tensor)是一个可用来表示在一些矢量,标量和其他张量之间的线性关系的多线性函数,这些线性关系的基本例子有内积,外积,线性映射以及笛卡儿积.其坐标在 n 维空间内,有$ n^r $ 个分量的一种量,其中每个分量都是坐标的函数,而在坐标变换时,这些分量也依照某些规则作线性变换.r称为该张量的秩或阶(与矩阵的秩和阶均无关系)

在同构的意义下,第零阶张量(r=0)即为标量,第一阶张量(r=1)即为矢量或者说向量,第二阶张量(r=2)则成为矩阵.由于变换方式的不同,张量分成"协变张量"(指标在下者),"逆变张量"(指标在上者),"混合张量"(指标在上和指标在下两者都有)三类.

一个张量的大小是一个向量,向量的第 n 个元素描述了张量在第 n 维上的大小,在numpy中,张量的大小即为shape

张量加减

张量加减依然是ufunc,这边不做复述.

Hadamard Product

Hadamard Product 要求形状一致,本质上也还是ufunc

A = np.array([
  [[1,2,3],    [4,5,6],    [7,8,9]],
  [[11,12,13], [14,15,16], [17,18,19]],
  [[21,22,23], [24,25,26], [27,28,29]],
  ])
B = np.array([
  [[1,2,3],    [4,5,6],    [7,8,9]],
  [[11,12,13], [14,15,16], [17,18,19]],
  [[21,22,23], [24,25,26], [27,28,29]],
  ])
C = A * B
C
array([[[  1,   4,   9],
        [ 16,  25,  36],
        [ 49,  64,  81]],

       [[121, 144, 169],
        [196, 225, 256],
        [289, 324, 361]],

       [[441, 484, 529],
        [576, 625, 676],
        [729, 784, 841]]])

张量的积算子

张量积算子通常表示为中间带有小x的圆$\otimes $.给定具有q维度的张量A和具有r维度的张量B，这些张量的乘积将是具有q + r的量级的新张量,或者换句话说q + r维度.

需要注意张量的积算子在r为1,2时是区别于向量点乘,矩阵点乘的.它要求形状一致

张量积算子不限于张量,也可以在矩阵和向量上执行，这可以是练习的好地方，以便发展更高维度的直觉.

在向量上:

$$ a = (a_1, a_2)$$

$$b = (b_1, b_2)$$

$$c = a \otimes b $$

$$c = {\begin{bmatrix} a_1 * b\\ a_2 * b \end{bmatrix} }$$

在矩阵上:

$$ A = {\begin{bmatrix} a_{1,1}&a_{1,2}\\ a_{2,1}&a_{2,2} \end{bmatrix}} $$

$$ B = {\begin{bmatrix} b_{1,1}& b_{1,2}\\ b_{2,1} & b_{2,2} \end{bmatrix}} $$

$$ C = A \otimes B $$

$$ C ={\begin{bmatrix} a_{1,1} * {\begin{bmatrix}b_{1,1}&b_{1,2}\b_{2,1}&b_{2,2}\end{bmatrix}}&a_{1,2}{\begin{bmatrix}b_{1,1}&b_{1,2}\b_{2,1}&b_{2,2}\end{bmatrix}}\ a_{2,1} * {\begin{bmatrix}b_{1,1}&b_{1,2}\b_{2,1}&b_{2,2}\end{bmatrix}}&a_{2,2}{\begin{bmatrix}b_{1,1}&b_{1,2}\b_{2,1}&b_{2,2}\end{bmatrix}} \end{bmatrix}} $$

$$ C = {\begin{bmatrix} a_{1,1} * b_{1,1}&a_{1,1} * b_{1,2}& a_{1,2} * b_{1,1}& a_{1,2} * b_{1,2}\\ a_{1,1} * b_{2,1}&a_{1,1} * b_{2,2}& a_{1,2} * b_{2,1}& a_{1,2} * b_{2,2}\\ a_{2,1} * b_{1,1}&a_{2,1} * b_{1,2}& a_{2,2} * b_{1,1}& a_{2,2} * b_{1,2}\\ a_{2,1} * b_{2,1}&a_{2,1} * b_{2,2}& a_{2,2} * b_{2,1}& a_{2,2} * b_{2,2} \end{bmatrix}} $$

注意上式中的值为压扁为矩阵后的形态,在拆解时每次拆解都是一层

在numpy中张量的积算子需要指定axes为0

np.tensordot(np.arange(4).reshape(2,2),np.arange(4).reshape(2,2),axes=0).shape
(2, 2, 2, 2)