フーリエ ニューラル演算子を使用した PDE の求解

R2026b
ライブ スクリプトを開く

この例では、偏微分方程式 (PDE) の解を出力するフーリエ ニューラル演算子 (FNO) ニューラル ネットワークに学習させる方法を示します。

"ニューラル演算子" [1] は、関数空間の間をマッピングするニューラル ネットワークの一種です。たとえば、システムの初期条件が与えられたときに PDE の解を出力するように学習できます。"フーリエ ニューラル演算子" (FNO) [2] は、フーリエ変換を活用することで、これらのマッピングをより効率的に学習できるニューラル演算子です。具体的には、空間領域における入力関数 f に対して、ニューラル ネットワークの層は、F(f)=fˆ となる変換 F を学習します。ここで fˆ は f のフーリエ変換です。これにより、FNO は周波数領域で動作できるようになります。ニューラル ネットワークを使用して PDE の解を予測すると、数値的な計算よりも高速に解を得られる場合があります。

次の図は、FNO を通るデータの流れを示しています。

Diagram of data flowing through FNO. The input is a function labeled "Initial Condition". The output is a function labeled "PDE Solution".

FNO は関数空間の間のマッピングを学習するため、学習データの解像度よりも高い解像度で出力関数を評価できます。この機能はゼロショット超解像 (ZSSR) として知られています。つまり、ZSSR は、追加の高解像度の例を必要とせずに、入力データと出力データの解像度を元の学習データを超えて向上させる機能です。次の図は、予測のための低解像度の学習データと高解像度の入力、およびそれらの出力の違いを示しています。

Plots showing low resolution and high resolution inputs and outputs. The low-resolution functions are composed of a small number of straight lines between points. The high-resolution functions are smooth curves.

この例では、1 次元バーガーズ方程式の解を出力するように FNO に学習させます。

学習データの読み込み

1 次元バーガーズ方程式のデータを読み込みます。

バーガーズ方程式は、流体力学、非線形音響学、気体力学、交通流など、応用数学のさまざまな分野で用いられる PDE です。その PDE は次のようなものです。

∂u∂t+u∂u∂x=1Re∂2u∂x2,x∈(0,1),t∈(0,1),

ここで、Re はレイノルズ数です。初期条件は次のとおりです。

u(x,0)=u0(x),x∈(0,1)

この方程式は、空間領域全体にわたって周期的境界条件を課します。目標は、初期状態 u0 から時刻 t=1 における解へのマッピング u0⟼u(x,1) を学習させることです。

1 次元バーガーズ方程式のデータを burgers_data_R10.mat から読み込みます。このデータには、バーガーズ方程式の初期速度 u0 と解 u(x,1) が含まれています。

dataFile = "burgers_data_R10.mat";
if ~exist(dataFile,"file")
    location = "https://ssd.mathworks.com/supportfiles/nnet/data/burgers1d/burgers_data_R10.mat";
    websave(dataFile,location); 
end

data = load(dataFile,"a","u");
u0 = data.a;
u1 = data.u;

学習データをいくつか表示します。

numPlots = 3;
gridSize = size(u0,2);
figure
tiledlayout(numPlots,2)
for i = 1:numPlots
    nexttile
    plot(linspace(0,1,gridSize),u0(i,:));
    title("Observation " + i + newline + "Initial Condition")
    xlabel("x")
    ylabel("u")

    nexttile
    plot(linspace(0,1,gridSize),u1(i,:));
    title("PDE Solution")
    xlabel("x")
    ylabel("v")
end

Three pairs of plots showing initial conditions (u vs. x) and corresponding PDE solutions (v vs. x) for observations 1, 2, and 3

学習データの準備

この例にサポート ファイルとして添付されている trainingPartitions 関数を使用して、学習データとテスト データを分割します。この関数にアクセスするには、例をライブ スクリプトとして開きます。データの 80% は学習に使用し、10% は検証に使用し、残りの 10% はテストに使用します。

numObservations = size(u0,1);
[idxTrain,idxValidation,idxTest] = trainingPartitions(numObservations,[0.8 0.1 0.1]);
U0Train = u0(idxTrain,:);
U1Train = u1(idxTrain,:);

U0Validation = u0(idxValidation,:);
U1Validation = u1(idxValidation,:);
U0Test = u0(idxTest,:);

FNO での処理量を減らすには、学習データと検証データをダウンサンプリングして、グリッド サイズが 1024 になるようにします。

gridSizeDownsampled = 1024;
downsampleFactor = floor(gridSize./gridSizeDownsampled);

U0Train = U0Train(:,1:downsampleFactor:end);
U1Train = U1Train(:,1:downsampleFactor:end);
U0Validation = U0Validation(:,1:downsampleFactor:end);
U1Validation = U1Validation(:,1:downsampleFactor:end);

FNO は、空間座標 (グリッド) に対応するチャネルをもつ入力データを必要とします。グリッドを作成し、入力データと連結します。

xMin = 0;
xMax = 1;
xGrid = linspace(xMin,xMax,gridSizeDownsampled);

xGridTrain = repmat(xGrid, [numel(idxTrain) 1]);
xGridValidation = repmat(xGrid, [numel(idxValidation) 1]);

U0Train = cat(3,U0Train,xGridTrain);
U0Validation = cat(3,U0Validation,xGridValidation);

学習セットと検証セットのサイズを表示します。1 番目、2 番目、および 3 番目の次元は、それぞれバッチ次元、空間次元、チャネル次元です。

size(U0Train)
ans = 1×3

    1638    1024    2


size(U0Validation)
ans = 1×3

    204    1024    2

ニューラル ネットワーク アーキテクチャの定義

fnoNetwork関数 (R2026b 以降) を使用して、FNO ネットワークを作成します。このネットワークは、リフティング層、直列接続された複数のフーリエ層、および投影層で構成されます。

Diagram of the neural network of layers connected in series. From left to right, there is the input, convolution, four fourier-ReLU blocks, convolution, ReLU, convolution, and output.

次のオプションを指定します。

  • 学習予測子 (初期状態) と目標値 (偏微分方程式の解) に合わせて入力サイズと出力サイズを指定します。

  • 1 次元 FNO ネットワークの場合、データが 1 つの空間次元をもつことを指定します。

  • 16 個のフーリエ モード、サイズが 64 である隠れ層、および ReLU 活性化をもつフーリエ層を 4 つ指定します。

inputSize = size(U0Train,3);
outputSize = size(U1Train,3);

net = fnoNetwork(inputSize,outputSize, ...
    Dimension=1, ...
    NumFourierLayers=4, ...
    NumModes=16, ...
    HiddenSize=64, ...
    Activation="relu");

R2026b より前: dlnetwork オブジェクトを使用して、ニューラル ネットワークのアーキテクチャを手動で定義します。

学習オプション

学習オプションを指定します。オプションの中から選択するには、経験的解析が必要です。実験を実行してさまざまな学習オプションの構成を調べるには、実験マネージャーアプリを使用できます。

  • SGDM ソルバーを使用して学習させます。

  • 初期学習率 0.001、低下係数 0.5 の区分的学習率スケジュールを使用して学習させます。

  • ミニバッチ サイズを 20 として 40 エポック学習させます。

  • すべてのエポックでデータをシャッフルします。

  • 入力データにはバッチ次元、空間次元、チャネル次元に対応する次元があるため、入力データ形式を "BSC" として指定します。

  • プロット内の学習の進行状況を監視します。

  • 検証データを指定し、100 回の反復ごとにネットワークを検証します。

  • 詳細出力を無効にします。

schedule = piecewiseLearnRate(DropFactor=0.5);

options = trainingOptions("sgdm", ...
    InitialLearnRate=1e-3, ...
    LearnRateSchedule=schedule, ...
    MaxEpochs=40, ...
    MiniBatchSize=20, ...
    Shuffle="every-epoch", ...
    InputDataFormats="BSC", ...
    Plots="training-progress", ...
    ValidationData={U0Validation,U1Validation}, ...
    ValidationFrequency=100, ...
    Verbose=false);

相対 L2 損失関数の定義

予測とターゲットを入力として受け取り、相対的な L2 損失を返す relativeL2Loss 関数を作成します。相対 L2 損失関数は、予測値とターゲットの差のユークリッド ノルムを使用してスケールを考慮し、それをターゲットのユークリッド ノルムで割る L2 損失関数のバリエーションです。

RelativeL2(Y,T)=∑n=1N‖Tn-Yn‖2‖Tn‖2

ここで、

  • Y は予測を表します。

  • T はターゲットを表します。

  • N は観測数を表します。

dlarray の入力から形式を削除するには、stripdims 関数を使用します。次元の順序を "SCB" として指定し、vecnorm で使用できるように空間次元を最初に配置します。

function loss = relativeL2Loss(Y,T)

Y = stripdims(Y,"SCB");
T = stripdims(T,"SCB");

p = vecnorm(T-Y,2,1);
q = vecnorm(T,2,1);

loss = sum(p./q);

end

ネットワークの学習

関数trainnetを使用してニューラル ネットワークに学習させます。カスタム損失関数 relativeL2Loss を使用します。既定では、trainnet 関数は利用可能な GPU がある場合にそれを使用します。GPU を使用するには、Parallel Computing Toolbox™ ライセンスとサポートされている GPU デバイスが必要です。サポートされているデバイスの詳細については、GPU 計算の要件 (Parallel Computing Toolbox)を参照してください。そうでない場合、関数は CPU を使用します。実行環境を手動で選択するには、ExecutionEnvironment 学習オプションを使用します。

net = trainnet(U0Train,U1Train,net,@relativeL2Loss,options);

Training progress plot showing training and validation loss decreasing over 3240 iterations across 40 epochs

ネットワークのテスト

真のラベルをもつホールドアウトされたテスト セットの予測を比較して、その平方根平均二乗誤差 (RMSE) を計算することで、モデルの予測パフォーマンスをテストします。

学習データの場合と同じ手順を使い、テスト データをダウンサンプリングして、グリッドを連結します。

xGridTest = repmat(xGrid, [numel(idxTest) 1]);

U0Test = U0Test(:,1:downsampleFactor:end);
U0Test = cat(3,U0Test,xGridTest);

U1Test = u1(idxTest,:);
U1Test = U1Test(:,1:downsampleFactor:end);

testnet関数を使用してニューラル ネットワークをテストします。既定では、testnet 関数は利用可能な GPU がある場合にそれを使用します。そうでない場合、関数は CPU を使用します。実行環境を手動で選択するには、testnet 関数の ExecutionEnvironment 引数を使用します。

rmseTest = testnet(net,U0Test,U1Test,"rmse",InputDataFormats="BSC")
rmseTest = 
0.0425

ZSSR の精度を評価するには、ダウンサンプリングせずにテスト データを使用してネットワークをテストします。

xGridZSSR = linspace(xMin,xMax,gridSize);
xGridTestZSSR = repmat(xGridZSSR, [numel(idxTest) 1]);

U0TestZSSR = u0(idxTest,:);
U0TestZSSR = cat(3,U0TestZSSR,xGridTestZSSR);
U1TestZSSR = u1(idxTest,:);

rmseTestZSSR = testnet(net,U0TestZSSR,U1TestZSSR,"rmse",InputDataFormats="BSC")
rmseTestZSSR = 
0.0425

新しいデータでの予測の実行

新しいデータを読み込みます。例として、最初の 3 つのテスト観測値を使用します。

dataNew = u0(idxTest(1:3),:);
numObservationsNew = size(dataNew,1);

グリッド サイズを計算し、それをデータと連結します。

gridSizeNew = size(dataNew,2);
xGridNew = linspace(xMin,xMax,gridSizeNew);
xGridNew = repmat(xGridNew, [size(dataNew,1) 1]);
U0New = cat(3,dataNew,xGridNew);

関数minibatchpredictを使用して予測を行います。既定では、minibatchpredict 関数は利用可能な GPU がある場合にそれを使用します。そうでない場合、関数は CPU を使用します。実行環境を手動で選択するには、ExecutionEnvironment オプションを使用します。

Y = minibatchpredict(net,U0New,InputDataFormats="BSC");

プロットで予測を可視化します。

figure
tiledlayout(numObservationsNew,2)
for i = 1:numObservationsNew
    nexttile
    plot(xGridNew(i,:),dataNew(i,:))
    title("Initial Condition")

    nexttile
    plot(xGridNew(i,:),Y(i,:))
    title("PDE Solution")
end

Three pairs of plots showing initial conditions and predicted PDE solutions for three test observations, with high-resolution output from the trained FNO

参考文献

  1. Li, Zongyi, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. “Neural Operator: Graph Kernel Network for Partial Differential Equations.” arXiv, March 6, 2020. http://arxiv.org/abs/2003.03485.

  2. Li, Zongyi, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. “Fourier Neural Operator for Parametric Partial Differential Equations,” 2020. https://openreview.net/forum?id=c8P9NQVtmnO.

参考

| | | | |

トピック