IDSS 复习笔记 | Computational Linear Algebra

基于自己的理解做的IDSS (Introduction to Data Science and Systems) 课程内容整理。 本文主要内容来自课程课件 Lecture 3: Complinear Algebra (v20202021b)。侧重点是代码实现和理论的结合。

邻接矩阵(Adjacency matrix)

    A B C D E F G H

A   0 1 1 1 0 0 0 0
B   0 0 0 1 1 0 0 0
C   0 1 0 1 0 0 0 0
D   0 0 0 0 0 1 1 1
E   0 0 0 0 0 0 0 0
F   0 0 0 0 0 0 0 0 
G   0 0 0 0 0 0 0 0
H   1 0 0 0 0 0 0 0
1
2
3
4
5
6
7
8
9
10
11
12
13
adj = np.array([[0, 1, 1, 1, 0, 0, 0, 0],
[0, 0, 0, 1, 1, 0, 0, 0],
[0, 1, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 1, 1],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 0, 0, 0, 0]])

in_degrees = np.sum(adj, axis=0)
out_degrees = np.sum(adj, axis=1)

np.allclose(adj, adj.T) #False #判断是否是无向图。
  • 顶点的入度(in-degree): 终点为该点的边数,是每一列的和。
  • 顶点的出度(out-degree): 起点为该点的边数,是每一行的和。
  • 无向图的邻接矩阵是对称的,即adj=adj.T
  • 将有向图变成无向图(使所有箭头变双向): \(A^\prime = A + A^T\)
  • 对于加权图,可以将图看做flows of "mass" through a network of vertices.
    • 如果一个顶点的总流出量\(>1\),即其行之和\(>1\),则称其是 a source of mass.
    • 如果一个顶点的总流出量\(<1\),即其行之和\(<1\),则称其是 a sink.
    • 如果一个顶点的总流出量就是\(1\),即其行之和为\(1\),then it conserves mass; it only ever re-routes things.

矩阵运算

矩阵的幂(Matrix powers/exponentiation)

\(C = A^n\),即\(n\)\(A\)相乘,\(A^4=AAAA\)

1
2
3
4
5
def powm(A, n):    
B = np.eye(A.shape[0])
for i in range(n):
B = A @ B
return B

稳定点(Stable point)

1
2
3
n_nodes = 30
adjacency = np.random.uniform(0, 1, (n_nodes, n_nodes))**50 #生成一个随机的变换矩阵
adjacency = (adjacency.T / np.sum(adjacency, axis=1)).T #按行归一化,使邻接矩阵是一个conserving matrix

1
2
3
4
5
6
7
8
9
10
11
12
13
timesteps = 60 #循环次数

test_vector = np.zeros((1, n_nodes)) #初始化一个行向量

node = np.random.randint(0, n_nodes) #设定一个随机的起始状态(将所有mass都放在某个“仓库”)
test_vector[0][node] = 1

vectors = [] #记录每一步后各个“仓库”的状态
vectors.append(test_vector.T) #按列向量存储

for i in range(timesteps):
test_vector = test_vector @ adjacency
vectors.append(test_vector.T)

  • 可以看到无论起始状态如何,最终都会达到稳定状态(steady state),而此时的向量正是邻接矩阵的特征向量之一,再次应用矩阵不会改变分布,即\({\bf x}_{t=59}={\bf x}_{t=60}\)

特征值(Eigenvalues)和特征向量(Eigenvectors)

\(A{\bf x} = \lambda {\bf x}\)

  • 特征向量: 矩阵进行线性变换时不会发生旋转,仅会缩放。
  • 特征值: 应用于特征向量的缩放比例。

幂迭代法(Power iteration method)求解

\({\bf x_n} = AAAA\dots AA{\bf x_0} = A^{n}{\bf x_0}\)

  • 直接迭代通常会导致数值爆炸或塌缩为零,需要在应用矩阵变换后对结果进行归一化来规避这一问题: \({\bf x_n} = \frac{A {\bf x_{n-1}}}{\|A{\bf x_{n-1}} \|_\infty}\)
  • 从效率考虑通常使用\(L_\infty\)进行归一化,也可以用别的范数。
  • 归一化的作用是消除缩放的影响,其他变换都不影响。
1
2
3
4
5
6
def power_iterate(A, x, n):        
for i in range(n):
x = A @ x #应用矩阵变换
x = x / np.linalg.norm(x, np.inf) #幂迭代时使用L-infinity范数进行归一化

return x / np.linalg.norm(x, 2) #最后除以L2范数,使可以与np.linalg.eig的结果进行比较(此条待完善)
1
2
3
4
5
6
A = np.random.normal(0, 1, (3, 3)) #生成一个随机的变换矩阵(3x3)
A = A + A.T #使其对称以保证每个特征值都是实数。->特征值为复数的邻接矩阵没有稳定状态,而是收敛为振荡状态。

random_vec = np.random.normal(0, 1, (3, 1)) #生成一个随机的初始向量

xn = power_iterate(A, random_vec, n=500) #对这个向量进行足够多次的幂运算以确保收敛到特征向量(leading eigenvector)

这一步求出的是前导特征向量(leading eigenvector),根据\(A{\bf x} = \lambda {\bf x}\) 可知特征值\(\lambda = \frac{A{\bf x}}{\bf x}\)

1
2
ratio = (A @ xn) / xn # elementwise ratio.shape=(3,1)
eigenvalue = ratio[0][0] # all elements of ratio will be the same

此时A @ xneigenvalue @ xn的值相等。

Numpy求解

1
evals, evecs = np.linalg.eig(A) #输出是特征值和特征向量
  • 一般来说,对于\(n \times n\)矩阵,np.linalg.eig将产生\(n\)个特征值和\(n\)列特征向量。
  • 特征向量是正交的(orthogonal),即任何一对特征向量的点积为零(此条待完善)
  • 对于非常大的矩阵,若只需求前导特征向量,幂迭代会比使用np.linalg.eig快很多(此条待完善)

一些其他相关

  • 特征向量仅存在于方阵中。
  • 对于任何矩阵,特征值都是唯一的,但特征向量不是唯一的(此条待完善)
  • np.linalg.eig由于float精度限制可能会导致数值不稳定,比如有时最小的特征向量不完全正交。
  • 如果矩阵满足某些特殊条件,则可以使用更稳定的算法。e.g.实数对称矩阵可以用np.linalg.eigh

本征谱(Eigenspectrum)

将特征值按大小排序的列表,特征值(的绝对值)大的一般认为更加“重要”。

1
evecs[:,np.argsort(-np.abs(evals))]

主成分分析(PCA,Principal Component Analysis)

此处需要复习一下均值、方差、协方差矩阵的概念。 协方差矩阵的特征向量就是主成分,指示了数据变化最大的方向。

1
2
3
4
5
6
x = np.random.normal(0, 1, (200, 2)) @ np.array([[0.05, 0.4], [-0.9, 1.0]]) #生成一个随机的数据集

evals, evecs = np.linalg.eig(np.cov(x, rowvar=False))

evals_sorted = evals[np.argsort(-np.abs(evals))]
evecs_sorted = evecs[:, np.argsort(-np.abs(evals))]

图中散点是随机生成的数据集,绿色和橙色是分离出来的两个主成分(即x的协方差的特征向量)椭圆是根据这两个值画的。三个圈只是为了好看,椭圆缩放比例是人工输入的。

特征分解(Eigendecomposition)

用特征向量和特征值来重构协方差矩阵: \(\Sigma = Q\Lambda Q^T\)。其中\(Q\)是单位特征向量(就是np.linalg.eig的输出),\(\Lambda\)是特征值的对角矩阵(对角线上是\(\lambda_i\),其余都为零)。

此处需要复习一下对角线矩阵的初始化方法。

1
2
3
4
xx = np.cov(x, rowvar=False)

evals, evecs = np.linalg.eig(np.cov(x, rowvar=False)) #evecs -> Q
reconstructed_A = evecs @ np.diag(evals) @ evecs.T

此时reconstructed_Axx的值相等。

  • 近似矩阵 处理高维数据时,\(\Sigma\)将是一个非常大的矩阵,我们希望可以仅存储前几个主成分,并通过这些数据来重构一个近似的\(\Sigma\)。近似的思路是截断特征值和特征向量(将小的特征值和特征向量设为0)。

    1
    2
    approx_evals = np.where(np.abs(evals)<1e-1, 0, evals) #将特征值中小于0.1的都截断为0
    approx_A = evecs @ np.diag(approx_evals) @ evecs.T #使用截断的特征值矩阵来重构A

  • 降维(Dimensionality reduction) 通过将原始数据集投影到几个协方差矩阵的主成分上(通常是较大的几个)对数据进行降维。 以2D数据为例,我们可以将原始的二维数据投影到较大的特征向量上以生成一维数据。(此条待完善)

    1
    x_projected = x @ evecs[:, 0] # remember x is a matrix with row vectors

矩阵的属性

矩阵的迹(Trace)

  • 方阵的迹就是对角线元素的和,即 \(\text{Tr}(A) = a_{1,1} + a_{2,2} + \dots + a_{n,n}\)
  • 矩阵的迹等于特征值的和: \(\text{Tr}(A) = \sum_{i=1}^n \lambda_i\)
  • 矩阵的迹可以看作是被应用矩阵做形变后parallelotope of a unit cube的周长。(严格来说是与周长成比例,比例系数是\(\text{Perimiter}(A)=2^{n-1} \text{Tr}(A)\))(此条待完善)

行列式(Determinant)

  • 行列式的值等于特征值的乘积: \(\text{det}(A) = \prod_{i=1}^n \lambda_i\)
  • 行列式可以看作是被应用矩阵做形变后parallelotope of a unit cube的体积,衡量线性变换后扩张或收缩的程度。(此条待完善)
  • 行列式值为零:
    • \(\text{det}(A)\)至少有一个特征值为零。
    • \(\text{det}(A)\)不可逆(因为在形变时会对向量空间进行一个或多个维度的折叠,导致信息丢失,所以该变换不可逆)。
    • \(\text{det}(A)\)奇异值矩阵

定矩阵(Definite matrix)和半定矩阵(Semi-definite matrix)

  • positive definite : 所有特征值都大于零,\(\lambda_i>0\)
  • positive semi-definite : 所有特征值都大于等于零,\(\lambda_i \geq0\)
  • negative definite : 所有特征值都小于零,\(\lambda_i<0\)
  • negative semi-definite : 所有特征值都小于等于零,\(\lambda_i \leq0\)
  • 正定矩阵(positive definite matrix)的属性: 对于任何非零的\({\bf x}\),正定矩阵\(A\)\({\bf x}^TA{\bf x}>0\)。(此条待完善)
    • 这表明\({\bf x}\)\(A {\bf x}\)的点积为正。
    • \({\bf x}\)\(A {\bf x}\)之间的夹角\(\theta < \frac{\pi}{2}\),因为\({\bf x}^TA{\bf x} = |{\bf x}||A {\bf x}|\cos{\theta}>0\) \(\rightarrow\) \(\cos{\theta}>0\)
    • 说明\(A\)\({\bf x}\)的旋转操作不超过 \(90^{\circ}\)
  • 正定性在讨论协方差矩阵和Hessian矩阵时有很重要的作用。

奇异值矩阵(Singular matrices)

\(\text{det}(A)=0\)的矩阵就是奇异值矩阵,不能求逆。可逆的矩阵是非奇异值矩阵。 如图,行列式为零则应用该矩阵时会使数据降维,导致信息丢失,所以不可逆。

1
2
3
#这是一个奇异值矩阵
sing_mat = np.array([[0.0, 1],
[0, 0.0]])

应用该矩阵进行变换后,二维数据降成了一维:

矩阵的秩(Rank)

矩阵的秩即非零奇异值的数量。

  • 满秩(full rank): 非零奇异值的数量=矩阵大小。
  • 满秩矩阵可逆且行列式不为零。
  • 不是满秩的矩阵是奇异值矩阵,不可逆且是deficient rank的。
  • 低秩(low rank): 非零奇异值的数量\(\ll\)矩阵大小

条件数(Condition number)

矩阵的条件数是最大奇异值和最小奇异值的比值。

  • 此定义仅针对满秩矩阵
  • 条件数衡量矩阵对微小变化的敏感度。
  • 条件数小的矩阵是well-conditioned的,不太会引起数值问题。
  • 条件数大的矩阵是ill-conditioned的,需要重视数值问题。
  • ill-conditioned几乎是奇异的,所以对其求逆可能会因为浮点数误差的原因导致无效。

矩阵的秩和条件数扩展了奇异的概念。

矩阵的秩可以看做是用来衡量矩阵“奇异程度”的值,在转换中损失了多少维度。条件数可以看做是衡量非奇异矩阵与奇异矩阵接近程度的值,接近奇异的矩阵可能会变得有效奇异。

矩阵的逆(Matrix inversion)

  • \(A^{-1}(A{\bf x}) = {\bf x}\),
  • \(A^{-1}A = I\)
  • \((A^{-1})^{-1} = A\)
  • \((AB)^{-1} = B^{-1}A^{-1}\)
  • 矩阵除法等效于左乘矩阵的逆。
  • 只有方阵可以求逆

计算矩阵的逆

1
2
A = np.random.normal(0, 1, (3, 3))
A_inv = np.linalg.inv(A)
1
2
3
x = np.random.normal(0, 1, (3, 1))
Ax = A @ x
x_ = np.linalg.inv(A) @ Ax

xx_的值相等。

特殊情况求逆

正交矩阵(orthogonal matrix)

  • 求逆: \(A^{-1}= A^T\)
  • 时间复杂度\(O(1)\)
  • 正交矩阵满足\(A^T=A^{-1}\)
  • 特征向量彼此正交(互相垂直,内积为零)
  • 特征值为1或-1
  • 任何纯旋转矩阵都是正交矩阵

对角矩阵(diagonal matrix)

  • 求逆: \(A^{-1} = \frac{1}{A}\)
  • 时间复杂度\(O(n)\)
  • 对角矩阵的所有非对角线元素为零
  • 对角矩阵的逆也是对角矩阵,其对角线元素是原始矩阵每个对角线元素的倒数
1
2
3
d = np.diag([-1, 2, -3, 4])
d_r = np.diag(1/np.diag(d))
d_inv = np.linalg.inv(d)

d_r=d_inv

稀疏矩阵(sparse matrices)

稀疏矩阵一般不求逆,因为它的逆通常不是稀疏矩阵,会消耗巨大的空间来储存。

伪逆(Pseudo-inverse)

当矩阵\(A\)不是正方形时,也可以通过求其伪逆\(A^+\)来得到近似可以“撤销”\(A\)变换的逆矩阵。伪逆可以利用SVD分解求逆来计算: \(A^{+} = V \Sigma^{-1} U^T\),或者使用Numpy的函数np.linalg.pinv

1
2
3
4
5
6
7
8
#创建一个非正方形的矩阵A
A = np.array([[0.5, 1.0, 2.0],
[1.0, 0.5, 0.0]])

y = np.array([[2], [5]])

A_pinv = np.linalg.pinv(A)# 求A的伪逆
x = A_pinv @ y # 用伪逆求解的x

  • A @ x\(\approx\)y
  • 非对称的情况下,np.linalg.pinv会找到能使\(\|A{\bf x} - {\bf y}\|_2\)最小化的\({\bf x}\)

线性方程组(Linear system of equations)

求解线性方程组

  • 对于\(A{\bf x} = {\bf y}\),如果\(A\)是正方形的,则两边同时左乘\(A^{-1}\)即可。这意味着\({\bf x}\)\({\bf y}\)需要是相同的形状。

  • 实践中通常不会用上面这个方法来求解线性系统,对高维矩阵求逆会使结果不稳定。
  • 常用的方法是优化迭代求解,通过最小化\(\|A{\bf x}-{\bf y}\|_2^2\)来得到特定的近似解(而不是求出所有解)。

奇异值分解(Singular value decomposition)

  • 特征分解仅适用于可对角化的方形矩阵(此条待完善)
  • 奇异值分解可以应用于不是方形的矩阵,更具有普适性。
  • \(A = U \Sigma V^T\)
    • \(A\)是任意矩阵,\(\mathbb{R}^{m \times n}\)
    • \(U\)方形的Unitary matrix(酉矩阵),\(\mathbb{R}^{m \times m}\),列向量包含left singular vectors
    • \(V\)是方形的Unitary matrix(酉矩阵),\(\mathbb{R}^{n \times n}\),列向量包含right singular vectors
      • Unitary matrix(酉矩阵): 共轭转置等于它的逆。\((\overline{U})^T=U^{-1}\)(此条待完善)
    • \(\Sigma\)是对角线矩阵,\(\mathbb{R}^{m \times n}\),对角线元素就是奇异值(Singular value)。
      • 奇异值与特征值密切相关但不完全相同,特殊情况除外。奇异值都是正实数。
  • \(A\)是实数矩阵时,\(U\)\(V\)都是正交的,且其行向量和列向量都是单位向量。
    • 此时\(U\)\(V\)可以看做纯旋转矩阵,而对角线矩阵\(\Sigma\)则是纯缩放矩阵。
    • SVD将任何矩阵变换都分解为旋转-缩放-旋转。
1
2
3
4
5
6
7
8
9
10
11
#创建一个非方形的矩阵A
A = np.array([[0.1, 0.2, 0.3],
[0.4, 0.5, 0.6]])

u, sigma, vt = np.linalg.svd(A)

#以下四个向量输出都是[1.0,1.0,1.0],U、V的行和列都是单位向量
np.linalg.norm(u, axis=1)
np.linalg.norm(u, axis=0)
np.linalg.norm(vt.T, axis=1)
np.linalg.norm(vt.T, axis=0)

u @ np.diag(sigma) @ vtA的值相等。

特殊情况: 对称半正定矩阵

  • 特征向量是\(U\)\(V\)的列向量(\(U\)\(V\)相同)
  • 特征值是\(\Sigma\)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
A = np.array([[2, 1, 0],
[1, 3, 1],
[0, 1, 4]])

evals, evecs = np.linalg.eig(A)

print("------------------")
print_matrix("\lambda (unordered)", evals)
print_matrix("\\bf{x} (unordered)", evecs)
print_matrix("\lambda", evals[::-1]) # Reverse order of eigenvalues to get them in size order
print_matrix("\\bf{x}", np.fliplr(evecs)) # Flip eigenvector matrix left right to correspond with ordered eigenvalues

u, sigma, vt = np.linalg.svd(A)
print_matrix("\Sigma", sigma)
print_matrix("U", u)
print_matrix("V", vt.T)
  • evals[::-1]sigma的值相等。(特征值逆序以使其按大小排序)
  • uvt.Tnp.fliplr(evecs)的值相等。(特征向量水平翻转以对应特征值)

用SVD求分数幂

\(A^n = V \Sigma^n U^T\) 比如求\(A^{1/2}\),就可以用\(A^{1/2} = U \Sigma^{1/2} V^T\)。(注意\(A_{1/2}\)不是单纯把\(A\)的每个元素开平方)

1
2
3
4
def powm(A, n): # generalised matrix power
u, sigma, vt = np.linalg.svd(A)
sigma_n = np.diag(sigma**n)
return u @ sigma_n @ vt

用SVD求逆

\(A ^{-1} = V \Sigma^{-1} U^T\)

1
2
3
4
5
6
7
8
9
10
11
def invert_by_svd(u, sigma, vt):
sigma_inv = np.zeros((u.shape[0], vt.shape[0]))
sqr = min(vt.shape[0], u.shape[0])
sigma_inv[:sqr, :sqr] = np.diag(1.0/sigma)
return vt.T @ sigma_inv @ u.T

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

u, sigma, vt = np.linalg.svd(A)

invert_by_svd(u, sigma, vt)np.linalg.inv(A)的输出相等。