簡介面向本碩博及科研人員這份Matlab仿真資源圍繞高斯過程與貝葉斯全局優(yōu)化方法展開包含完整可運行的算法代碼與演示錄像適合用于貝葉斯優(yōu)化編程學(xué)習(xí)、課程設(shè)計、畢業(yè)論文驗證及課題組算法預(yù)研。資源共11個文件以9個Matlab腳本.m為核心另有1個文本格式的說明文檔和1個操作演示視頻整體壓縮包僅173KB輕量緊湊、便于按需查看。目前已有2031人學(xué)習(xí)下載屬于高頻使用的教研資料。內(nèi)容上代碼模塊化實現(xiàn)了高斯過程回歸、期望改進(jìn)EI、置信上界UCB等關(guān)鍵步驟綜合運行腳本可一鍵串聯(lián)優(yōu)化流程說明文檔補充了必要的算法注釋與運行提示操作視頻則直觀展示從主腳本運行、子函數(shù)調(diào)用到當(dāng)前路徑設(shè)置的完整細(xì)節(jié)可有效規(guī)避直接運行子函數(shù)等常見錯誤幫助讀者快速掌握貝葉斯全局優(yōu)化在Matlab中的落地寫法。1. 貝葉斯全局優(yōu)化把調(diào)參黑匣子變成一條可復(fù)現(xiàn)的收斂曲線貝葉斯全局優(yōu)化Bayesian OptimizationBO在 MATLAB 里最常見的打開方式是bayesopt一行命令但大多數(shù)教程止步于“填個函數(shù)、按回車、看曲線”它背后的高斯過程代理模型、采集函數(shù)選點、新樣本回填這三步到底怎么協(xié)作很多人沒機會看細(xì)節(jié)。這次拆的資源包包含一套手動實現(xiàn)的 BO 仿真代碼和配套操作視頻能把黑匣子一層層剝開每一步怎么算、每個參數(shù)調(diào)了有什么影響、為什么幾十次迭代就能逼近全局最優(yōu)。適合正在做超參數(shù)搜索、仿真試驗次數(shù)有限、或者想從“會用bayesopt”進(jìn)階到“能自己改 BO 邏輯”的 MATLAB 用戶。下面的內(nèi)容按原理、復(fù)現(xiàn)、避坑、驗證往下走代碼可以直接抄。2. 高斯過程代理模型與采集函數(shù)BO 少算的那幾步憑什么是“全局最優(yōu)”2.1 高斯過程回歸在 BO 里的角色均值、協(xié)方差與超參數(shù)BO 的核心假設(shè)是目標(biāo)函數(shù)很貴我不能隨便算。所以先用一個代理模型模仿它代理模型給每個未知點同時輸出“預(yù)測值”和“不確定度”。工程上最常用的代理模型就是高斯過程回歸Gaussian Process RegressionGPR。區(qū)別在于普通回歸只給一個y_predGP 給的是y_pred ± 若干倍標(biāo)準(zhǔn)差這個標(biāo)準(zhǔn)差會隨訓(xùn)練點密度的變化而變化離已知樣本越遠(yuǎn)標(biāo)準(zhǔn)差越大也就是“我心里越?jīng)]底”。在 MATLAB 里fitrgp就是干這件事的。它有三個東西需要你理解均值函數(shù)默認(rèn)是常數(shù)也就是整體趨勢。對 BO 場景常數(shù)均值夠用因為采集函數(shù)選點主要靠協(xié)方差帶來的空間結(jié)構(gòu)。協(xié)方差函數(shù)核函數(shù)決定“兩個點距離多近才算相似”。平方指數(shù)核平滑但偏“軟”Matern 5/2 核在工程目標(biāo)函數(shù)上更穩(wěn)因為它對參數(shù)變化不那么敏感不容易把響應(yīng)面擬合得過于震蕩。噪聲參數(shù)Sigma目標(biāo)函數(shù)如果是仿真軟件每次結(jié)果可能帶數(shù)值噪聲如果是確定性數(shù)學(xué)函數(shù)Sigma可以設(shè)得很小。參數(shù)調(diào)優(yōu)時fitrgp默認(rèn)用最大似然估計去擬合超參數(shù)包括核函數(shù)的長度尺度、噪聲方差。這意味著每一輪 BO 迭代其實都在重新學(xué)一次“響應(yīng)面的形狀”。這也是 BO 耗時的大頭后面避坑章會細(xì)說。核函數(shù)MATLAB 中名稱適用目標(biāo)函數(shù)特點常用場景平方指數(shù)squaredexponential非常光滑、無突變理論測試函數(shù)Matern 5/2matern52存在輕微粗糙、有噪聲工程仿真、CFD/ FEMMatern 3/2matern32波動較強結(jié)構(gòu)優(yōu)化響應(yīng)面ARD 類ardsquaredexponential不同維度敏感度差異大維度多、需自動降權(quán)我一般做 BO 仿真第一版都用 Matern 5/2不是因為它理論最漂亮而是實際跑工程目標(biāo)函數(shù)時它比平方指數(shù)更不容易出現(xiàn)那種“過度自信”的窄尖峰假象。這一條在后面避坑章會反復(fù)出現(xiàn)。2.2 采集函數(shù) EI怎么平衡“挖已知最優(yōu)點”和“探索未知區(qū)”代理模型建好了下一步要決定下一個試驗點在哪里。這個決策函數(shù)叫采集函數(shù)Acquisition Function。最經(jīng)典的是期望改進(jìn)Expected ImprovementEI已知當(dāng)前已觀測到的最小值是f_best某候選點x的 GP 預(yù)測均值為mu標(biāo)準(zhǔn)差為sigma該點的期望改進(jìn)量約為EI (f_best - mu) * Phi(z) sigma * phi(z)其中z (f_best - mu) / sigma直觀理解如果mu比f_best小很多說明這個點很有潛力如果sigma很大說明這里沒被探索過也可能藏著更優(yōu)解。EI 同時考慮這兩項所以它既不會死守當(dāng)前最優(yōu)也不會漫無目的地亂跑。這就是 BO“全局”二字的來源。采集函數(shù)不只有 EI 一種。UCB 是把mu和sigma做加權(quán)和偏好更激進(jìn)地探索POI 則只關(guān)心“超過當(dāng)前最優(yōu)的概率”容易保守。MATLAB 的bayesopt里可以直接指定AcquisitionFunctionName手動實現(xiàn)時我通常選 EI因為它參數(shù)少、對噪聲容忍度更高。如果你想更保守可以把 EI 公式里的f_best換成f_best xixi取 0.01 到 0.1這樣會給探索加一點權(quán)重避免卡在局部極值。2.3 BO 與網(wǎng)格搜索、隨機搜索的邊界在哪方法采樣次數(shù)對低維參數(shù)效果好對高維參數(shù)效果好能否利用歷史信息網(wǎng)格搜索指數(shù)增長是否否隨機搜索固定次數(shù)否是否BO固定次數(shù)是弱是網(wǎng)格搜索在 2 個參數(shù)時還能忍到 6 個參數(shù)時每維度取 10 個點就是 10^6 次仿真顯然不現(xiàn)實。隨機搜索雖然在高維下比網(wǎng)格好用但每一輪都是獨立抽簽之前算出的好點不會指導(dǎo)后面往哪個方向補采。BO 的價值就是“每一輪都站在過去所有采樣點肩膀上”。但 BO 不是萬能。參數(shù)維度超過 20 且沒有明顯低維結(jié)構(gòu)時GP 的協(xié)方差矩陣訓(xùn)練成本和采集函數(shù)優(yōu)化成本都會陡增這時 BO 反而不如隨機搜索加后處理。這個邊界要心里有數(shù)BO 適合的是“單次評估貴、預(yù)算小、維度在 2 到 15 之間”的優(yōu)化問題正好覆蓋大多數(shù)仿真調(diào)參場景。3. 在 MATLAB 里搭建 BO 仿真從測試函數(shù)到完整迭代循環(huán)3.1 用 Branin-Hoo 標(biāo)準(zhǔn)測試函數(shù)當(dāng)目標(biāo)為什么選它為了驗證實現(xiàn)正確需要選一個全局最優(yōu)已知、形狀又不太簡單的目標(biāo)函數(shù)。Branin-Hoo 是全局優(yōu)化領(lǐng)域的標(biāo)準(zhǔn)測試函數(shù)定義域x1 ∈ [-5, 10]、x2 ∈ [0, 15]有三個全局最優(yōu)點最小值約為0.3979。它有兩個局部極小能考驗 BO 是否會過早收斂。工程里的仿真黑箱通常也有多個局部極值因此這個測試函數(shù)非常有代表性。下面是初始采樣與目標(biāo)函數(shù)定義% Branin-Hoo 目標(biāo)函數(shù)x 是 1x2 行向量 obj (x) (x(2) - 5.1./(4*pi^2).*x(1).^2 5/pi.*x(1) - 6).^2 ... 10*(1 - 1./(8*pi)).*cos(x(1)) 10; lb [-5, 0]; % 兩個維度的下界 ub [10, 15]; % 兩個維度的上界 % 用 Latin Hypercube 生成 5 個初始樣本保證空間不扎堆 rng(42); X0 lb lhsdesign(5, 2) .* (ub - lb); y0 zeros(5, 1); for i 1:5 y0(i) obj(X0(i, :)); end這里lhsdesign生成的是[0,1]區(qū)間內(nèi)的拉丁超立方樣本通過.*(ub - lb) lb縮放到實際定義域。初始樣本數(shù)量不要太多5 個就夠因為 BO 的價值就是在迭代中逐步補點。rng(42)固定隨機種子保證復(fù)現(xiàn)時初始點完全一致。再看 GP 擬合與迭代循環(huán)X X0; y y0; for iter 1:30 % 1. 用當(dāng)前所有觀測點擬合高斯過程 gpr fitrgp(X, y, ... KernelFunction, matern52, ... Sigma, 1e-3, ... SigmaLowerBound, 1e-5, ... Standardize, true); % 2. 在候選網(wǎng)格上預(yù)測均值和標(biāo)準(zhǔn)差 n 100; x1g linspace(lb(1), ub(1), n); x2g linspace(lb(2), ub(2), n); [X1, X2] ndgrid(x1g, x2g); Xcand [X1(:), X2(:)]; [mu, sigma] predict(gpr, Xcand); sigma max(sigma, 1e-6); % 防止除零 % 3. 計算 EI 并挑出最大期望改進(jìn)點 fbest min(y); z (fbest - mu) ./ sigma; EI (fbest - mu) .* normcdf(z) sigma .* normpdf(z); [~, idx] max(EI); xnew Xcand(idx, :); ynew obj(xnew); % 4. 把新點加入歷史樣本進(jìn)入下一輪 X [X; xnew]; y [y; ynew]; end這段代碼的邏輯很直白每次迭代先重新用fitrgp擬合 GP然后在 100×100 的網(wǎng)格候選點上算 EI取 EI 最大的那個點作為本次仿真點算完真實目標(biāo)值后回填到樣本庫。fitrgp的Sigma參數(shù)是高斯過程噪聲標(biāo)準(zhǔn)差SigmaLowerBound給一個下限防止fitrgp把噪聲調(diào)整為負(fù)數(shù)或過小導(dǎo)致數(shù)值奇異。Standardize設(shè)為true會對輸出做標(biāo)準(zhǔn)化對 Branin 這種取值范圍跨度大的目標(biāo)函數(shù)很有必要。normcdf和normpdf分別是標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)和概率密度函數(shù)對應(yīng)前文 EI 公式里的Phi(z)和phi(z)。網(wǎng)格密度n100在二維問題下生成一萬個候選點計算量完全可接受。如果你優(yōu)化的是三維建議降到 40否則候選笛卡爾積會爆炸到 6 萬以上。3.2 用 MATLAB 內(nèi)置 bayesopt 走一遍結(jié)果對照用手動實現(xiàn)適合理解原理、改采集函數(shù)、加自定義核函數(shù)。但在實際工程里你更可能直接用內(nèi)置bayesopt它內(nèi)部做了更多數(shù)值優(yōu)化比如采集函數(shù)的最優(yōu)解不是靠網(wǎng)格找的而是用連續(xù)優(yōu)化器去細(xì)搜因此結(jié)果會更精細(xì)。% 定義可優(yōu)化變量及其邊界 optVar [optimizableVariable(x1, [-5, 10]); ... optimizableVariable(x2, [0, 15])]; % 直接調(diào)用內(nèi)置貝葉斯優(yōu)化 results bayesopt(obj, optVar, ... AcquisitionFunctionName, expected-improvement, ... MaxObjectiveEvaluations, 30, ... IsObjectiveDeterministic, true, ... Verbose, 1); [Xopt, yopt] bestPoint(results, Criterion, min-observed); fprintf(最優(yōu)值: %.4f 在 x1%.2f x2%.2f\n, yopt, Xopt.x1, Xopt.x2);MaxObjectiveEvaluations控制總仿真次數(shù)包括初始點。IsObjectiveDeterministic設(shè)為true告訴優(yōu)化器函數(shù)無噪聲這樣它會直接插值而不額外估計噪聲如果你的目標(biāo)函數(shù)每次運行有一點點浮點誤差或仿真隨機性要改成false。Verbose設(shè)為 1 能在命令行實時看到每輪迭代結(jié)果。跑完以后bestPoint的min-observed是找“已觀測到的最小值”比min-expected更穩(wěn)妥后者是用代理模型預(yù)測的可能出現(xiàn)偏差。把手動循環(huán)和bayesopt的結(jié)果放在一起對比你會發(fā)現(xiàn)最優(yōu)值都很接近0.3979但bayesopt通常用更少的迭代數(shù)找到更接近最優(yōu)的位置因為它在選點上做了連續(xù)優(yōu)化而不是粗網(wǎng)格。不過手動版的價值恰恰在“可改”你想換成 UCB 采集函數(shù)、加一個禁止采樣的禁區(qū)或者把 GP 換成帶梯度信息的變體手動版才是起點。3.3 收斂曲線怎么畫、怎么看BO 的收斂曲線不應(yīng)該畫“歷史最優(yōu)值”因為最小值會一直單調(diào)下降顯得很漂亮但信息量不夠。我建議畫兩個子圖% 子圖1當(dāng)前最優(yōu)值隨評估次數(shù)變化 figure; plot(1:length(y), min(y), b-, LineWidth, 1.5); xlabel(評估次數(shù)); ylabel(當(dāng)前最優(yōu)目標(biāo)值); % 子圖2每次新增樣本的改進(jìn)量 dy [0; min(y(2:end)) - min(y(1:end-1))]; hold on; bar(1:length(dy), dy); legend(當(dāng)前最優(yōu), 單次改進(jìn)量);單次改進(jìn)量如果長期小于1e-4說明 BO 已經(jīng)收斂或卡住了。這時要看最后兩個樣本點是否落在一起如果落在同一個局部極值周圍且 EI 最大值依然很高說明探索不足需要調(diào)大xi或改用 UCB。3.4 資源包里操作視頻對應(yīng)的復(fù)現(xiàn)節(jié)奏這套資源里附帶了一段操作視頻。視頻里按“生成初始樣本 → 跑第一輪 GP 擬合 → 觀察 EI 面 → 加入新樣本 → 再看收斂曲線”的順序演示對應(yīng)上面代碼的執(zhí)行路徑。錄制操作時建議打開 MATLAB 的編輯器與命令行分區(qū)每輪迭代之間停頓 3 到 5 秒讓EI面的變化過程在視頻里可讀。你拿到手復(fù)現(xiàn)時只要固定rng(42)每一步的輸出數(shù)字會和視頻完全一致可以用來核對你的 MATLAB 版本和工具箱是否裝齊。注意fitrgp屬于 Statistics and Machine Learning Toolboxbayesopt屬于 Global Optimization Toolbox。缺少任一個工具箱都會在調(diào)用時報未定義函數(shù)或許可證錯誤。4. BO 仿真避坑核矩陣奇異、邊界徘徊與超參翻車的排查記錄4.1 現(xiàn)象擬合 GP 時 MATLAB 報錯“無法計算協(xié)方差矩陣”運行fitrgp時直接報錯有時還提示Matrix is singular或者covariance matrix not positive definite。新手最容易以為是fitrgp不穩(wěn)定其實是樣本里有重復(fù)點或者核函數(shù)長度尺度被訓(xùn)練得極小導(dǎo)致兩點之間協(xié)方差趨近零整個矩陣接近奇異。原因BO 迭代時如果候選網(wǎng)格太粗EI 最大值點可能會反復(fù)選到同一位置或者初始樣本點之間距離太近造成數(shù)值重復(fù)。遇到這類問題首先要做去重檢查而不是急著換核函數(shù)。解決在進(jìn)入下一輪前判斷新點與已有樣本的最小距離小于閾值就跳過本次迭代。或者在fitrgp里給Sigma加一個下限比如1e-5這相當(dāng)于給協(xié)方差矩陣對角線加噪聲數(shù)值穩(wěn)定很多。實際操作中我會在循環(huán)里加一行% 檢查重復(fù)點 dmin min(sqrt(sum((X - xnew).^2, 2))); if dmin 1e-6 continue; end4.2 現(xiàn)象EI 最大值一直貼著邊界新樣本全在參數(shù)上下限處這可能是 BO 最隱蔽的翻車現(xiàn)場曲線也在下降但選出來的點永遠(yuǎn)是x2的上界或下界。第一次遇到會懷疑采集函數(shù)寫錯了其實多半是目標(biāo)函數(shù)在邊界外的趨勢沒被 GP 學(xué)準(zhǔn)或者網(wǎng)格候選點沒有把邊界附近加密。原因Branin 這類函數(shù)的最優(yōu)點在邊界內(nèi)側(cè)但如果 GP 長度尺度被最大似然估計得偏大邊界外會被預(yù)測成同樣低的值EI 就會把點推到邊界上。另一個常見原因是Standardize沒開GP 對超出訓(xùn)練范圍的外推值過于自信。解決先看訓(xùn)練好的核函數(shù)長度尺度如果兩個維度的長度尺度都遠(yuǎn)遠(yuǎn)大于定義域?qū)挾日f明模型把函數(shù)學(xué)成了近似線性需要給長度尺度設(shè)置上限。手動實現(xiàn)時我用的是KernelFunction, matern52并且把候選網(wǎng)格在邊界處加密比如linspace改用兩側(cè)加密的分布。更實際的辦法是改用bayesopt它會用連續(xù)優(yōu)化采集函數(shù)而不受網(wǎng)格分辨率限制邊界問題會輕很多。4.3 現(xiàn)象每次迭代都要幾十秒30 次根本跑不完fitrgp每次迭代都在重新估計全部超參數(shù)數(shù)據(jù)點從 5 漲到 35雖然樣本量不大但每次訓(xùn)練都要多次計算核矩陣的逆和似然函數(shù)梯度累積起來就很慢。如果你把候選網(wǎng)格設(shè)成n300一次預(yù)測就是 9 萬個點更慢。原因手動實現(xiàn)里把“GP 擬合”放在了每一輪循環(huán)內(nèi)而實際上超參數(shù)在樣本量變化不大時沒必要每輪都重學(xué)。解決一種常見做法是每 5 輪更新一次超參數(shù)其余輪次用上一輪訓(xùn)練好的gpr對象直接預(yù)測。MATLAB 的predict可以用舊模型不強制重新擬合。把fitrgp調(diào)用放進(jìn)if mod(iter, 5) 1條件里能省一半以上的訓(xùn)練時間。同時把候選網(wǎng)格降到n80二維就是 6400 個點對選點精度影響很小。4.4 現(xiàn)象目標(biāo)函數(shù)帶噪聲時BO 結(jié)果反而比隨機搜索差工程仿真里目標(biāo)函數(shù)常有數(shù)值噪聲比如有限元結(jié)果的小幅波動、CFD 迭代未完全收斂。直接用手動版跑EI 會把噪聲當(dāng)成真信號導(dǎo)致它反復(fù)去采樣那些“碰巧低”的點。原因EI 公式里的sigma是 GP 對模型不確定度的估計而fitrgp默認(rèn)會嘗試把噪聲Sigma擬合到接近下限模型只會去解釋觀測值的波動而不是平滑噪聲。目標(biāo)函數(shù)有噪聲時這會推動 BO 進(jìn)入過度擬合狀態(tài)。解決fitrgp里把Sigma初值設(shè)成目標(biāo)值方差的 10%并設(shè)置SigmaLowerBound為一個明顯不為零的值比如sqrt(var(y)) * 0.01。如果噪聲隨時間均勻還可以在同一參數(shù)點重復(fù)采樣三次取均值。最省事的路徑是直接用bayesopt并把IsObjectiveDeterministic設(shè)為false它會自動估計一個合適的噪聲方差。4.5 現(xiàn)象照視頻復(fù)現(xiàn)收斂曲線對不上資源包的操作視頻里每一步都有命令行輸出但自己跑出來的最優(yōu)值和曲線走勢不一致。這時八成不是代碼抄錯而是隨機種子、版本差異或初始采樣差異。原因視頻里用rng(42)固定了初始點如果你跳過了這一行l(wèi)hsdesign每次生成的 5 個點都不一樣。另一個常見因素是高版本 MATLAB 的fitrgp默認(rèn)迭代次數(shù)或優(yōu)化容差變了導(dǎo)致同一批數(shù)據(jù)擬合出的超參數(shù)略有差異。解決在腳本最前面寫清rng(42)復(fù)現(xiàn)時把fitrgp的OptimizeHyperparameters默認(rèn)行為關(guān)掉比較麻煩。更通用的是把每輪X、y保存成.mat文件視頻里也做同樣記錄對比中間過程數(shù)據(jù)而不是只對比最終圖。save(sprintf(bo_iter_%02d.mat, iter), X, y);5. 用已知最優(yōu)反推實現(xiàn)正確性采樣預(yù)算壓縮與結(jié)果驗證拿到這套代碼第一件事不是立刻換自己的工程目標(biāo)函數(shù)而是先在 Branin-Hoo 上確認(rèn)實現(xiàn)正確。定義域內(nèi)已知全局最優(yōu)值約為0.3979三組最優(yōu)解分別落在(-pi, 12.275)、(pi, 2.275)、(9.42478, 2.475)附近。如果能跑到0.398左右說明代理模型和 EI 實現(xiàn)沒問題如果卡在某個局部最優(yōu)附近比如0.5以上且 30 輪不再下降就要回頭查采集函數(shù)里的符號是不是寫反了。我會做一個更嚴(yán)格的驗證把MaxObjectiveEvaluations分別壓到 10、20、30跑三次各自記錄最終最優(yōu)值和第一次低于0.42的輪次。結(jié)果應(yīng)該呈現(xiàn)一個規(guī)律——10 次通常能到0.45左右20 次能穩(wěn)定低于0.4130 次逼近0.398。如果 10 次反而比 30 次結(jié)果更好那說明迭代中出現(xiàn)了數(shù)值抖動或重復(fù)采樣需要檢查去重邏輯。驗證完正確性后再把這個 BO 封裝成自己的函數(shù)輸入是目標(biāo)函數(shù)句柄、邊界、迭代次數(shù)輸出是全局最優(yōu)解、最優(yōu)值和歷史樣本記錄。用 MATLAB 的 OOP 架構(gòu)把測試函數(shù)、采集函數(shù)和代理模型拆成三個類后續(xù)換到自己工程目標(biāo)時就不用來回改腳本了。我現(xiàn)在每換一個新仿真場景都會先固定隨機種子在 2 維標(biāo)準(zhǔn)函數(shù)上跑通一輪再套真實仿真目標(biāo)這一條習(xí)慣救過我很多次。希望幫到你。本文還有配套的精品資源點擊獲取