沒有採樣的高維推理問題中的不確定性估計？
我正在研究一個高維推理問題（大約 2000 個模型參數），我們能夠通過結合使用基於梯度的優化和遺傳算法找到對數後驗的全局最大值來穩健地執行 MAP 估計。
除了找到 MAP 估計值之外，我非常希望能夠對模型參數的不確定性做出一些估計。
我們能夠有效地計算對數後驗相對於參數的梯度，所以從長遠來看，我們的目標是使用哈密頓量 MCMC 進行一些採樣，但現在我對基於非採樣的估計感興趣。
我知道的唯一方法是在模式下計算 Hessian 的逆，以將後驗近似為多元正態，但即使這樣對於如此大的系統似乎也不可行，因為即使我們計算 $ \sim 4\times10^{6} $ 黑森州的元素我敢肯定我們找不到它的逆。
誰能建議在這種情況下通常使用哪種方法？
謝謝！
編輯- 有關問題的其他信息
背景
這是一個與大型物理實驗相關的逆問題。我們有一個描述一些物理場的二維三角形網格，我們的模型參數是網格每個頂點處這些場的物理值。網格有大約 650 個頂點，我們對 3 個場進行建模，這就是我們的 2000 個模型參數的來源。
我們的實驗數據來自不直接測量這些場的儀器，而是這些場的複雜非線性函數的量。對於每種不同的儀器，我們都有一個前向模型，它將模型參數映射到實驗數據的預測，並且預測和測量之間的比較產生對數似然。
然後，我們總結了所有這些不同工具的對數似然，並添加了一些對數先驗值，這些值對這些字段應用了一些物理約束。
因此，我懷疑這個“模型”是否完全屬於一個類別——我們無法選擇模型是什麼，它取決於收集我們實驗數據的實際儀器的功能。
數據集
數據集由 500x500 幅圖像組成，每個攝像頭有一張圖像，因此總數據點為 500x500x4 = $ 10^6 $ .
誤差模型
我們現在將問題中的所有誤差都視為高斯誤差。在某些時候，我可能會嘗試轉移到 student-t 誤差模型，只是為了獲得一些額外的靈活性，但僅使用 Gaussians 似乎仍然可以正常工作。
似然示例
這是一個等離子體物理實驗，我們的絕大多數數據來自指向等離子體的相機，鏡頭前有特定的過濾器，僅查看光譜的特定部分。
重現數據有兩個步驟；首先，我們必須對來自網格上的等離子體的光進行建模，然後我們必須將該光建模回相機圖像。
不幸的是，對來自等離子體的光進行建模取決於什麼是有效的速率係數，即在給定場的情況下，不同過程會發出多少光。這些速率是由一些昂貴的數值模型預測的，因此我們必須將它們的輸出存儲在網格上，然後進行插值以查找值。速率函數數據只計算一次——我們存儲它，然後在代碼啟動時從它構建一個樣條，然後該樣條被用於所有函數評估。
認為 $ R_1 $ 和 $ R_2 $ 是速率函數（我們通過插值評估），然後是發射 $ i $ ‘網格的第一個頂點 $ \mathcal{E}_i $ 是（誰）給的 $$ \mathcal{E}_i = R_1(x_i, y_i) + z_i R_2(x_i, y_i) $$ 在哪裡 $ (x,y,z) $ 是我們在網格上建模的 3 個場。將發射向量獲取到相機圖像很容易，只需與矩陣相乘 $ \mathbf{G} $ 它編碼了每個相機像素所看到的網格的哪些部分。
由於誤差是高斯的，所以這個特定相機的對數似然是 $$ \mathcal{L} = -\frac{1}{2} (\mathbf{G}\vec{\mathcal{E}} - \vec{d})^{\top}\mathbf{\Sigma}^{-1} (\mathbf{G}\vec{\mathcal{E}} - \vec{d}) $$
在哪裡 $ \vec{d} $ 是相機數據。總對數似然是上述 4 個表達式的總和，但對於不同的相機，它們都有不同版本的速率函數 $ R_1, R_2 $ 因為他們正在觀察光譜的不同部分。
先驗示例
我們有各種先驗，它們實際上只是為各種數量設置了某些上限和下限，但這些先驗往往不會對問題產生太大影響。我們確實有一個作用很強的先驗，它有效地將拉普拉斯類型平滑應用於場。它也採用高斯形式： $$ \text{log-prior} = -\frac{1}{2}\vec{x}^{\top}\mathbf{S}\vec{x} -\frac{1}{2}\vec{y}^{\top}\mathbf{S}\vec{y} -\frac{1}{2}\vec{z}^{\top}\mathbf{S}\vec{z} $$
首先，我認為您的統計模型是錯誤的。我將您的符號更改為統計學家更熟悉的符號，因此讓
$$ \mathbf{d}=\mathbf{y}=(y_1,\dots,y_N),\ N=10^6 $$
成為您的觀察（數據）向量，並且
$$ \begin{align} \mathbf{x}&=\boldsymbol{\theta}=(\theta_1,\dots,\theta_p) \ \mathbf{y}&=\boldsymbol{\phi}=(\phi_1,\dots,\phi_p) \ \mathbf{z}&=\boldsymbol{\rho}=(\rho_1,\dots,\rho_p), \ p \approx 650 \ \end{align} $$
你的參數向量，總維度 $ d=3p \approx 2000 $ . 然後，如果我理解正確，你假設一個模型
$$ \mathbf{y} = \mathbf{G}\mathbf{r_1}(\boldsymbol{\theta}, \boldsymbol{\phi})+\boldsymbol{\rho}\mathbf{G}\mathbf{r_2}(\boldsymbol{\theta}, \boldsymbol{\phi}))+\boldsymbol{\epsilon},\ \boldsymbol{\epsilon}\sim\mathcal{N}(0,I_N) $$
在哪裡 $ \mathbf{G} $ 是個 $ N\times d $ 樣條插值矩陣。
這顯然是錯誤的。來自同一相機的圖像中不同點處的誤差，以及來自不同相機的圖像中相同點處的誤差不可能是獨立的。您應該研究空間統計和模型，例如廣義最小二乘法、半變異函數估計、克里金法、高斯過程等。
話雖如此，由於您的問題不是模型是否是實際數據生成過程的良好近似，而是如何估計這樣的模型，我將向您展示一些選項來做到這一點。
HMC
2000 個參數並不是一個很大的模型，除非你是在筆記本電腦上訓練這個東西。數據集更大（ $ 10^6 $ 數據點），但是，如果您可以訪問云實例或帶有 GPU 的機器，那麼Pyro或Tensorflow Probability等框架將很快解決此類問題。因此，您可以簡單地使用 GPU 驅動的 Hamiltonian Monte Carlo。
優點：“精確”推斷，在鏈中無限數量的樣本的限制下。
缺點：對估計誤差沒有嚴格限制，存在多個收斂診斷指標，但沒有一個是理想的。
大樣本逼近
濫用符號，讓我們表示 $ \theta $ 通過連接三個參數向量獲得的向量。然後，使用貝葉斯中心極限定理 (Bernstein-von Mises)，您可以近似 $ p(\theta\vert \mathbf{y}) $ 和 $ \mathcal{N}(\hat{\theta_0}_n,I_n^{-1}(\theta_0)) $ ， 在哪裡 $ \theta_0 $ 是“真”參數值， $ \hat{\theta_0}_n $ 是 MLE 估計 $ \theta_0 $ 和 $ I_n^{-1}(\theta_0) $ 是費雪信息矩陣評估在 $ \theta_0 $ . 當然， $ \theta_0 $ 未知，我們將使用 $ I_n^{-1}(\hat{\theta_0}_n) $ 反而。Bernstein-von Mises 定理的有效性取決於您可以在這裡找到的一些假設，例如：在您的情況下，假設 $ R_1,R_2 $ 是光滑和可微的，這個定理是有效的，因為高斯先驗的支持是整個參數空間。或者，更好的是，如果您的數據實際上是您所假設的 iid，那將是有效的，但我不相信它們是，正如我在開始時解釋的那樣。
優點：在 $ p<<N $ 案子。保證收斂到正確的答案，在獨立同分佈的情況下，當似然是平滑和可微的並且先驗在鄰域中是非零的 $ \theta_0 $ .
缺點：正如您所指出的，最大的缺點是需要反轉 Fisher 信息矩陣。此外，我不知道如何憑經驗判斷近似值的準確性，除非使用 MCMC 採樣器從 $ p(\theta\vert \mathbf{y}) $ . 當然，這首先會破壞使用 B-vM 的實用性。
變分推理
在這種情況下，而不是找到確切的 $ p(\theta\vert \mathbf{y}) $ （這需要計算一個 $ d- $ i 維積分），我們選擇近似 $ p $ 和 $ q_{\phi}(\theta) $ ， 在哪裡 $ q $ 屬於參數族 $ \mathcal{Q}_{\phi} $ 由參數向量索引 $ \phi $ . 我們尋找 $ \phi^* $ st 之間的一些差異度量 $ q $ 和 $ p $ 被最小化。選擇這個度量作為 KL 散度，我們得到變分推理方法：
$$ \DeclareMathOperator*{\argmin}{arg,min} \phi^*=\argmin_{\phi\in\Phi}D_{KL}(q_{\phi}(\theta)||p(\theta\vert\mathbf{y})) $$
要求 $ q_{\phi}(\theta) $ :
- 它應該是可區分的 $ \phi $ ，以便我們可以應用隨機梯度下降等大規模優化方法來解決最小化問題。
- 它應該足夠靈活，可以準確地近似 $ p(\theta\vert\mathbf{y}) $ 對於一些價值 $ \phi $ ，但也足夠簡單，很容易從中取樣。這是因為估計 KL 散度（我們的優化目標）需要估計期望值 $ q $ .
你可能會選擇 $ q_{\phi}(\theta) $ 完全因式分解，即 $ d $ 單變量概率分佈：
$$ q_{\phi}(\theta)=\prod_{i=1}^d q_{\phi_i}(\theta_i) $$
這就是所謂的平均場變分貝葉斯方法。可以證明（例如，參見本書的第 10 章）每個因素的最優解 $ q_{\phi_j}(\theta_j) $ 是
$$ \log{q_j^*(\theta_j)} = \mathbb{E}_{i\neq j}[\log{p(\mathbf{y},\theta)}] + \text{const.} $$
在哪裡 $ p(\mathbf{y},\theta) $ 是參數和數據的聯合分佈（在您的情況下，它是您的高斯似然和參數上的高斯先驗的乘積）並且期望是相對於其他變分單變量分佈 $ q_1^(\theta_1),\dots,q_{j-1}^(\theta_{j-1}),q_{j+1}^(\theta_{j+1}),\dots,q_{d}^(\theta_{d}) $ . 當然，由於其中一個因素的解決方案取決於所有其他因素，我們必須應用迭代過程，初始化所有分佈 $ q_{i}(\theta_{i}) $ 進行一些初始猜測，然後使用上面的等式一次迭代地更新它們。請注意，不是將上面的期望計算為 $ (d-1)- $ 維積分，這在先驗和似然不共軛的情況下會令人望而卻步，您可以使用蒙特卡羅估計來近似期望。
平均場變分貝葉斯算法不是您可以使用的唯一可能的 VI 算法：Kingma & Welling, 2014 中提出的變分自動編碼器，“自動編碼變分貝葉斯”是一個有趣的替代方案，其中，而不是假設完全分解的形式為了 $ q $ ，然後為 $ q_i $ , $ q $ 被假定為多元高斯，但在每個可能具有不同的參數 $ N $ 數據點。為了分攤推理成本，使用神經網絡將輸入空間映射到變分參數空間。有關算法的詳細描述，請參閱論文：VAE 實現再次在所有主要的深度學習框架中可用。