标签：count 重构 bar 数据 模型 线性 cp hat
什么是网络重构
本文省略大量复杂网络、图重构预备知识，并以极简化的方式介绍基于网络博弈产生的动力学数据进行线性图重构的过程。
网络重构的意义
大数据中有一类数据是由复杂网络代表的实际动力学系统产生的。
由拓扑结构中各部分生成的这部分数据可以被测量，但对产生此类数据的结构我们一无所知。
所以进行网络结构的预测可以让我们更好的对复杂系统进行理解和控制。
不过网络重构绝非易事，有四个基本难点需要克服：
- 拓扑复杂性
- 复杂系统的高度非线性
- 动力学数据存在不可避免的噪声
- 收集数据时有效信息可能部分缺失
-
拓扑复杂性不难理解，图中任意节点均会受到其他 k-hop 邻居的影响。且有向图和无向图代表的信息不相同；同质性和异质性可能影响重构精度；以及其他拓扑结构的特殊性，均会导致重构不准。
-
实际系统的高度非线性正是建模的最大难点，由于图中的信息传递是多路径的、同时的、局部或全局的，所以在系统上进行的动力学过程输出的序列数据将会是高度序列化和无法先验感知的，也会是高度 “Nonlinear” 的。
-
噪声就更好理解了，任何以预测作为目的的统计方法均带有误差和偏差。噪声伴随收集到的数据“与生俱来”，任意重构模型的误差大体分为外生误差、内生误差，其中外生误差可以通过某些手段加以控制，而内生误差几乎无法避免，所以建模时需要格外谨慎处理。
-
最后是数据的缺失问题。在我们收集到的数据中，总会存在数据缺失或无效的情况，而后来的一些研究采用各种奇特而精妙的手段解决了这个问题，这种情况不包含在我们所讨论的线性重构模型中。
这种网络重构不止停留在理论研究中，同样有落地的实例，比如最近的两个医学方面的有趣例子：
- 学者 Marko Jusup 利用临床数据重构大脑认知过程和调控机理。
www.youtube.com/watch?v=lew6ApzPQsg
- 基于医学数据的 MRI 影像重建，本质上是一种线性重构模型，后来又引入机器学习相关方法。
Hammernik K, Klatzer T, Kobler E, et al. Learning a variational network for reconstruction of accelerated MRI data[J]. Magnetic resonance in medicine, 2018, 79(6): 3055-3071.
线性重构的动力学原理
现实中很多的复杂动力学过程都可以用一个简单的微分方程表示为：
\[\dfrac{dx(t)}{dt} = F(x(t)) + \Tau (t) \]其中 \(x(t) = (x_1(t), x_2(t), \cdots, x_N(t))\) 表示拓扑结构中 N 个节点的状态变量。
\(F(x(t)) = (F_1(x(t)), F_2(x(t)), \cdots, F_N(x(t)))\) 表示除节点本事之外的高阶动力学过程。
\(\Tau(t) = (\Tau_1(t), \Tau_2(t), \cdots, \Tau_N(t))\) 用于刻画动力学过程受噪声的影响。
然而由于现实中刻画复杂性并不容易，找到一个最优的 \(F(x(t))\) 十分困难，但可以通过最简单的线性形式近似代替它，即：
\[\dfrac{dx(t)}{dt} = \hat{A}x(t) + \Tau (t) \]其中，定义待估邻接矩阵 A 为 \([A_{ij}]\), 此时通过观测得到的全部节点输出数据 \(x(t) = (x_1(t), x_2(t), \cdots, x_N(t))\) ,可以近似表示线性表达式与邻接矩阵之间的关联。关联矩阵定义为：
\[\hat{C} = [C_{ij}], C_{ij} = <x_i(t)x_j(t)> = \dfrac{1}{T} \int_0^T x_i(t) x_j(t) dt \]其中T取足够长时间的数据，这样通过 Fokker-Planck 方程的定义，可以将上述微分方程简化为：
\[\hat{A}\hat{C} + \hat{C}\hat{A}^{T} = -\hat{Q} \]其中 \(\hat{Q} = [Q_{ij}]\) 是一个常数矩阵。该式的问题是已知 x(t) 的时间序列后不能直接求出常数矩阵 A 和 Q 的值。
关于 Fokker-Planck 方程定义与求解：Risken H. The Fokker-Planck Equation. Heidelberg: Springer-Verlag, 1984.
而且实际问题中 \(x_i(t)\) 不可能连续测量, 所以只能假设可测得离散的数据：\(x_i(t_1),x_i(t_2),\cdots, x_i(t_L), i=1,2,\cdots,N,L>1\)
假设以固定频率测量变量数据，则：
\[\Delta t = t_{k+1} - t_k,k=1,\cdots,L-1 \]满足 \(0 <t_{k+1} - t_k <1\)
此时可以给出测量数据变量变化的速率：
\[\dot{x}_i(t_k) = \dfrac{x_i(t_{k+1})-x_i(t_k)}{\Delta t},k=1,\cdots,L-1 \]于是可以利用测量数据和加性白噪声假设求出邻接矩阵 A 和常数矩阵 Q 的值（线性近似表达式）：
\[\dot{x}(t_k) = \hat{A}x(t_k)+\Tau(t^{'}_k) \]对这个等式进行进一步变形，两遍右乘 \(x(t)^{T}\):
\[\hat{B} = \hat{A}\hat{C}+<\Tau (t')x(t)^T> \]其中：\(\hat{B} = [B_{ij}]\), \(B_{ij} = \dfrac{1}{L-1}\sum_{k=1}^{L-1}\dot{x}_i(t_k)x_j(t_k)\)
离散情况下 \(C_{ij} = \dfrac{1}{L-1}\sum_{k=1}^{L-1} x_i(t_k)x_j(t_k)\)
在线性重构时一般假设加性白噪声和节点动力学数据在统计上不相关，所以有：
\[<\Tau (t')x(t)^T> = 0 \]如此一来，就有了动力学重构的基本线性求解表达式：
\[\hat{B} = \hat{A}\hat{C} \]其中，B一般为收集到的动力学现象数据，C通常为与B相关联的过程数据，也可以被观测到。A为待估的邻接矩阵，也就是我们最终要的网络结构。
模型求解
对上述的线性模型求解的根本目的是估计邻接矩阵\hat{A}，该矩阵是一个二元矩阵，即元素非零即一。所以如何在求解时将估计参数有效压缩到0和1成为求解时首要考虑的问题。
早期对这种线性模型一般采取压缩感知（Compress sensing）的方法，即 Penalty function 为：
\[\mathop{min}\limits_{A}||B-AC||_1 ,s.t. \quad B=AC \]但这种线性感知的方法无法很好的对估计范围进行压缩，所以自然的引入惩罚项，让压缩范围朝 0-1 靠近。于是最早在2015年的 PRL 上发表了一种改进的 Lasso 方法，它的Penalty function为：
\[\mathop{min}\limits_{A}\left[\dfrac{1}{2}||B-AC||_2^2 + \lambda ||A||_1 \right],s.t. \quad B=AC \]最后回忆一下动力学数据生成，也就是重构中的 B 和 C，见下图：
基于Lasso方法的线性重构求解
原创，仅供自己回忆使用，请勿转载
此代码作为一个求解模块，配合随笔【理解与实现基于网络博弈的动力学数据生成】中的代码使用。
For Python
# Lasso solve
import cvxpy as cp
X_bar = []
for i in range(G.number_of_nodes()):
X = cp.Variable(G.number_of_nodes())
# obj = cp.Minimize(cp.multiply((1 / (2 * L)), cp.square(cp.norm2((np.array(Y[:,0]).reshape(L, )) - (cp.matmul(Phi[i], X))))) + cp.multiply(.001 ,cp.norm1(X)))
obj = cp.Minimize(cp.multiply((1 / (2 * L)), cp.square(cp.norm2(Y_test[:, i] - (cp.matmul(Phi[i], X))))) + cp.multiply(.001, cp.norm1(X)))
prob = cp.Problem(obj)
prob.solve(solver = cp.CVXOPT) # ['CVXOPT', 'ECOS', 'ECOS_BB', 'GLPK', 'GLPK_MI', 'OSQP', 'SCS']
X_bar.append(X.value)
# print("\n\n status: {} \n\n optimal value: {} \n\n optimal varible: x = {}".format(prob.status, prob.value, X.value))
X_bar = np.matrix(X_bar) # not be exactly zero or one
# test
count_nl, count_el, count_fail = 0, 0, 0
for i in range(X_bar.shape[0]):
for j in range(X_bar.shape[1]):
if X_bar[i, j] > -.1 and X_bar[i, j] < .1:
count_nl += 1
# X_bar[i, j] = 0
elif X_bar[i, j] > .9 and X_bar[i, j] < 1.1:
count_el += 1
# X_bar[i, j] = 1
else: count_fail += 1
pre_fail.append(count_fail / (count_fail + count_el + count_nl))
np.savetxt('A_M_bar.csv', X_bar, delimiter=',')
end2 = time()
A_M = np.loadtxt(open("C:\\Users\\27225\\Desktop\\workshop\\A_M.csv", "rb"), delimiter=",", skiprows=0)
A_M_bar = np.loadtxt(open("C:\\Users\\27225\\Desktop\\workshop\\A_M_bar.csv", "rb"), delimiter=",", skiprows=0)
# a = np.array((A_M.flatten(), A_M_bar.flatten()))
# a = a.T[np.lexsort(-a)].T
# compute the AUROC
fpr, tpr, _ = roc_curve(A_M.flatten(), A_M_bar.flatten())
roc_auc = auc(fpr, tpr)
print("AUROC = {}".format(roc_auc))
plt.plot(fpr, tpr, color = 'darkorange', lw = 2, label = 'ROC curve (AUROC = %0.2f)' % roc_auc)
plt.legend(loc = 'best')
标签：count,重构,bar,数据,模型,线性,cp,hat 来源： https://www.cnblogs.com/huaiyu-DL-dynamics/p/15982243.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。