【網路重構】理解與實現基於線性模型的圖重構
什麼是網路重構
本文省略大量複雜網路、圖重構預備知識，並以極簡化的方式介紹基於網路博弈產生的動力學資料進行線性圖重構的過程。
網路重構的意義
大資料中有一類資料是由複雜網路代表的實際動力學系統產生的。
由拓撲結構中各部分生成的這部分資料可以被測量，但對產生此類資料的結構我們一無所知。
所以進行網路結構的預測可以讓我們更好的對複雜系統進行理解和控制。
不過網路重構絕非易事，有四個基本難點需要克服：
- 拓撲復雜性
- 複雜系統的高度非線性
- 動力學資料存在不可避免的噪聲
- 收集資料時有效資訊可能部分缺失
-
拓撲復雜性不難理解，圖中任意節點均會受到其他 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')