📜用例
📜C++和Python通信引文道路社评电商大规模行为图结构数据模型 | 📜应用于MNIST、瑞士卷、双峰数据集 | 📜算法:多维缩放、等距映射、局部线性嵌入 、谱嵌入 | 📜耦合拉普拉斯算子生成对齐几何图形特征空间 | 📜分数拉普拉斯算子对非平滑回归函数的非参数回归方法处理 | 📜图拉普拉斯算子处理谱聚类
✒️Python稀疏矩阵离散拉普拉斯算子
在本文中,我们考虑使用 python/numpy/scipy 计算一维数组的离散拉普拉斯算子。我们引入稀疏矩阵技术,该技术对于计算某些类型的扩散偏微分方程相当有效。
我们假设反应堆只包含四个单元,因此我们可以将物质向量写成以下形式:
v
⃗
=
(
v
1
v
2
v
3
v
4
)
\vec{v}=\left(\begin{array}{l} v_1 \\ v_2 \\ v_3 \\ v_4 \end{array}\right)
v=
v1v2v3v4
让我们将相邻像元之间的距离表示为 h 并使用周期性边界条件。然后离散拉普拉斯算子计算如下:
Δ
v
⃗
=
1
h
(
v
2
−
v
1
h
−
v
1
−
v
4
h
v
3
−
v
2
h
−
v
2
−
v
1
h
v
4
−
v
3
h
−
v
3
−
v
2
h
v
1
−
v
4
h
−
v
4
−
v
3
h
)
=
1
h
2
(
v
2
−
2
v
1
+
v
4
v
3
−
2
v
2
+
v
1
v
4
−
2
v
3
+
v
2
v
1
−
2
v
4
−
v
3
)
\Delta \vec{v}=\frac{1}{h}\left(\begin{array}{l} \frac{v_2-v_1}{h}-\frac{v_1-v_4}{h} \\ \frac{v_3-v_2}{h}-\frac{v_2-v_1}{h} \\ \frac{v_4-v_3}{h}-\frac{v_3-v_2}{h} \\ \frac{v_1-v_4}{h}-\frac{v_4-v_3}{h} \end{array}\right)=\frac{1}{h^2}\left(\begin{array}{l} v_2-2 v_1+v_4 \\ v_3-2 v_2+v_1 \\ v_4-2 v_3+v_2 \\ v_1-2 v_4-v_3 \end{array}\right)
Δv=h1
hv2−v1−hv1−v4hv3−v2−hv2−v1hv4−v3−hv3−v2hv1−v4−hv4−v3
=h21
v2−2v1+v4v3−2v2+v1v4−2v3+v2v1−2v4−v3
因此,拉普拉斯算子对离散空间中向量的作用相当于其乘以离散拉普拉斯矩阵:
Δ
v
⃗
=
1
h
2
(
−
2
1
0
1
1
−
2
1
0
0
1
−
2
1
1
0
1
−
2
)
v
⃗
\Delta \vec{v}=\frac{1}{h^2}\left(\begin{array}{cccc} -2 & 1 & 0 & 1 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ 1 & 0 & 1 & -2 \end{array}\right) \vec{v}
Δv=h21
−21011−21001−21101−2
v
离散一维拉普拉斯矩阵可以针对诺依曼边界条件进行简化:
Δ
=
1
h
2
(
−
1
1
0
0
1
−
2
1
0
0
1
−
2
1
0
0
1
−
1
)
\Delta=\frac{1}{h^2}\left(\begin{array}{cccc} -1 & 1 & 0 & 0 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ 0 & 0 & 1 & -1 \end{array}\right)
Δ=h21
−11001−21001−21001−1
或对于狄利克雷边界条件:
显然,可以轻松地为任意数量的单元构建这样的矩阵。
由于大量的零元素,离散拉普拉斯矩阵始终是稀疏的。因此我们可以使用稀疏矩阵技术。它在 scipy.sparse 库中实现,因此我们需要导入它:
import scipy.sparse as sp
稀疏矩阵可以以不同的格式存储。 您可以轻松地将一种转换为另一种。例如:
mm = sp.lil_matrix((10000,10000), dtype=float32)
mm[100,2345] = 12.5
通常稀疏矩阵的所有对角线元素都不为零(这正是我们的情况)。 在这种情况下,我们可以使用特殊函数 spdiags。 我们按如下方式对矩阵对角线进行编号:主对角线的编号为“0”,主对角线下方的对角线的编号为“-1”,主对角线上方的对角线的编号为“1”,依此类推。所有对角线都被认为是给定的 1xN 形式,但“-1”对角线向左移动,“1”向右移动。 例如,代码
diag1 = [1]*10000
diag2 = [2]*10000
mm = sp.spdiags([ diag2, diag1, diag2], [ -1, 0, 1], 10000, 10000)
形成一个三对角矩阵,主对角线为“1”,相邻对角线为“2”。接下来,我们可以显式设置一些单元格的值:
mm = mm.tolil()
mm[9999,0]=1
mm[0,9999]=1
当我们想要进一步执行向量矩阵乘法时,我们应该将稀疏矩阵转换为 CSR 形式:
mm = mm.tocsr()
m = np.random.randn(10000)
n = mm.dot(m)
让我们实现构建具有诺伊曼边界条件的离散拉普拉斯稀疏矩阵的代码:
def gen_laplacian(num):
data1 = [-2]*num
data1[0] = -1
data1[-1] = -1
data2 = [1]*num
L = sp.spdiags([ data2, data1, data2 ], [-1, 0, 1], num, num)
return L.tocsr()
主对角线由“-2”元素组成,我们将“-1”元素放在末尾。要应用周期性边界条件,我们应该将“1”包含在整个矩阵的左下角和右上角。我们可以使用 LIL 形式来做到这一点:
def gen_laplacian_pbc(num):
data1 = [-2]*num
data2 = [1]*num
L = sp.spdiags([ data2, data1, data2 ], [-1, 0, 1], num, num).tolil()
L[num-1,0] = 1
L[0,num-1] = 1
return L.tocsr()
现在可以应用拉普拉斯滤波器:
u = numpy.randn(10000)
L = gen_laplacian_pbc(10000)
u_transformed = L.dot(u)
或逆拉普拉斯滤波器:
u_inv = lg.spsolve(L,u)
请注意,要使用 spsolve,应导入稀疏矩阵的线性代数包。