L-BFGS 法と深層学習を使用した偏微分方程式の求解

この例では次を使用します。

この例では、記憶制限 BFGS (L-BFGS) アルゴリズムを使用してバーガース方程式の解を数値的に計算し、Physics Informed Neural Network (PINN) の学習を行う方法を説明します。

バーガース方程式は、応用数学のさまざまな分野に登場する偏微分方程式 (PDE) です。特に、流体力学、非線形音響学、気体力学、および交通流の分野です。L-BFGS アルゴリズム [1] は、Broyden-Fletcher-Goldfarb-Shanno (BFGS) アルゴリズムを近似する準ニュートン法です。

与えられた計算領域 $[- 1, 1] \times [0, 1]$ について、この例では、Physics Informed Neural Network (PINN) [2] を使用して、標本 $(x, t)$ を入力として受け取る多層パーセプトロンニューラルネットワークに学習させます。ここで、 $x \in [- 1, 1]$ は空間変数、 $t \in [0, 1]$ は時間変数です。また、 $u (x, t)$ (u は次のバーガース方程式の解) が返されます。

$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} - \frac{0.01}{π} \frac{\partial^{2} u}{\partial x^{2}} = 0,$

初期状態は $u (x, t = 0) = - s i n (π x)$ 、境界条件は $u (x = - 1, t) = 0$ および $u (x = 1, t) = 0$ とします。

この例では、与えられた入力 $(x, t)$ について、ネットワークの出力 $u (x, t)$ がバーガース方程式、境界条件、および初期状態を満たすように強制することで、モデルに学習させます。このモデルの学習には、事前のデータ収集を必要としません。PDE の定義および制約を使用して、データを生成できます。

学習データの生成

モデルに学習させるには、境界条件を強制し、初期状態を強制し、バーガース方程式を満たす選点のデータセットが必要です。

25 個の等間隔の時間点を選択して、各境界条件 $u (x = - 1, t) = 0$ および $u (x = 1, t) = 0$ を強制します。

numBoundaryConditionPoints = [25 25];

x0BC1 = -1*ones(1,numBoundaryConditionPoints(1));
x0BC2 = ones(1,numBoundaryConditionPoints(2));

t0BC1 = linspace(0,1,numBoundaryConditionPoints(1));
t0BC2 = linspace(0,1,numBoundaryConditionPoints(2));

u0BC1 = zeros(1,numBoundaryConditionPoints(1));
u0BC2 = zeros(1,numBoundaryConditionPoints(2));

50 個の等間隔の時間点を選択して、初期状態 $u (x, t = 0) = - s i n (π x)$ を強制します。

numInitialConditionPoints  = 50;

x0IC = linspace(-1,1,numInitialConditionPoints);
t0IC = zeros(1,numInitialConditionPoints);
u0IC = -sin(pi*x0IC);

初期状態および境界条件について、データをグループ化します。

X0 = [x0IC x0BC1 x0BC2];
T0 = [t0IC t0BC1 t0BC2];
U0 = [u0IC u0BC1 u0BC2];

10,000 個の点を選択し、ネットワークの出力がバーガース方程式を満たすように強制します。

numInternalCollocationPoints = 10000;

points = rand(numInternalCollocationPoints,2);

dataX = 2*points(:,1)-1;
dataT = points(:,2);

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

それぞれ 20 個の隠れニューロンがある 9 つの全結合演算で多層パーセプトロンニューラルネットワークアーキテクチャを定義します。最初の全結合演算には、入力 $x$ および $t$ に対応する 2 つの入力チャネルがあります。最後の全結合演算には、1 つの出力 $u (x, t)$ があります。

ネットワークのハイパーパラメーターを指定します。

numLayers = 9;
numNeurons = 20;

ニューラルネットワークを作成します。

layers = featureInputLayer(2);

for i = 1:numLayers-1
    layers = [
        layers
        fullyConnectedLayer(numNeurons)
        tanhLayer];
end

layers = [
    layers
    fullyConnectedLayer(1)]

layers = 
  18×1 Layer array with layers:

     1   ''   Feature Input     2 features
     2   ''   Fully Connected   20 fully connected layer
     3   ''   Tanh              Hyperbolic tangent
     4   ''   Fully Connected   20 fully connected layer
     5   ''   Tanh              Hyperbolic tangent
     6   ''   Fully Connected   20 fully connected layer
     7   ''   Tanh              Hyperbolic tangent
     8   ''   Fully Connected   20 fully connected layer
     9   ''   Tanh              Hyperbolic tangent
    10   ''   Fully Connected   20 fully connected layer
    11   ''   Tanh              Hyperbolic tangent
    12   ''   Fully Connected   20 fully connected layer
    13   ''   Tanh              Hyperbolic tangent
    14   ''   Fully Connected   20 fully connected layer
    15   ''   Tanh              Hyperbolic tangent
    16   ''   Fully Connected   20 fully connected layer
    17   ''   Tanh              Hyperbolic tangent
    18   ''   Fully Connected   1 fully connected layer

層配列を dlnetwork オブジェクトに変換します。

net = dlnetwork(layers)

net = 
  dlnetwork with properties:

         Layers: [18×1 nnet.cnn.layer.Layer]
    Connections: [17×2 table]
     Learnables: [18×3 table]
          State: [0×3 table]
     InputNames: {'input'}
    OutputNames: {'fc_9'}
    Initialized: 1

  View summary with summary.

学習可能パラメーターのデータ型が double である場合、PINN の学習の精度が向上します。関数 dlupdate を使用して、ネットワークの学習可能パラメーターを double に変換します。たとえば single 型の学習可能パラメーターに依存する GPU 最適化を使用するネットワークなど、すべてのニューラルネットワークが double 型の学習可能パラメーターをサポートしているわけではないことに注意してください。

net = dlupdate(@double,net);

モデル損失関数の定義

この例の最後のモデル損失関数のセクションにリストされている関数 modelLoss を作成します。この関数は、ニューラルネットワーク、ネットワーク入力、初期状態、および境界条件を入力として受け取り、学習可能なパラメーターについての損失の勾配、および対応する損失を返します。

学習オプションの指定

学習オプションを指定します。

1500 回反復して学習させる
勾配またはステップのノルムが 0.00001 より小さい場合は、学習を早期に停止する。
L-BFGS ソルバーの状態に対して既定のオプションを使用する。

maxIterations = 1500;
gradientTolerance = 1e-5;
stepTolerance = 1e-5;
solverState = lbfgsState;

ニューラルネットワークの学習

学習データを dlarray オブジェクトに変換します。入力 X および T の形式を "BC" (batch、channel) と指定し、初期状態の形式を "CB" (channel、batch) と指定します。

X = dlarray(dataX,"BC");
T = dlarray(dataT,"BC");
X0 = dlarray(X0,"CB");
T0 = dlarray(T0,"CB");
U0 = dlarray(U0,"CB");

dlaccelerate を使用して、損失関数を高速化します。

accfun = dlaccelerate(@modelLoss);

L-BFGS の更新に関する損失関数を含む関数ハンドルを作成します。自動微分を使用して関数 modelLoss 内で関数 dlgradient を評価するには、関数 dlfeval を使用します。

lossFcn = @(net) dlfeval(accfun,net,X,T,X0,T0,U0);

TrainingProgressMonitor オブジェクトを初期化します。各反復で損失をプロットし、勾配とステップのノルムを監視します。monitor オブジェクトを作成するとタイマーが開始されるため、学習ループに近いところでオブジェクトを作成するようにしてください。

monitor = trainingProgressMonitor( ...
    Metrics="TrainingLoss", ...
    Info=["Iteration" "GradientsNorm" "StepNorm"], ...
    XLabel="Iteration");

各反復でデータセット全体を使用し、カスタム学習ループを使用してネットワークに学習させます。それぞれの反復で次を行います。

関数 lbfgsupdate を使用して、ネットワークの学習可能なパラメーターとソルバーの状態を更新します。
学習の進行状況モニターを更新します。
勾配のノルムが指定された許容値を下回る場合、または直線探索が失敗した場合は、学習を早期に停止します。

iteration = 0;
while iteration < maxIterations && ~monitor.Stop
    iteration = iteration + 1;
    [net, solverState] = lbfgsupdate(net,lossFcn,solverState);

    updateInfo(monitor, ...
        Iteration=iteration, ...
        GradientsNorm=solverState.GradientsNorm, ...
        StepNorm=solverState.StepNorm);

    recordMetrics(monitor,iteration,TrainingLoss=solverState.Loss);

    monitor.Progress = 100 * iteration/maxIterations;

    if solverState.GradientsNorm < gradientTolerance || ...
            solverState.StepNorm < stepTolerance || ...
            solverState.LineSearchStatus == "failed"
        break
    end

end

モデルの精度の評価

$t$ の値が 0.25、0.5、0.75、および 1 である場合について、相対 $l^{2}$ 誤差を使用して、深層学習モデルの予測値とバーガース方程式の真の解を比較します。

モデルをテストする際のターゲット時間を設定します。各時間について、[-1,1] の範囲内の 1001 個の等間隔の点で解を計算します。

tTest = [0.25 0.5 0.75 1];
numPredictions = 1001;
XTest = linspace(-1,1,numPredictions);
XTest = dlarray(XTest,"CB");

モデルをテストします。

figure
tiledlayout("flow")

for i=1:numel(tTest)
    t = tTest(i);
    TTest = t*ones(1,numPredictions);
    TTest = dlarray(TTest,"CB");

    % Make predictions.
    XTTest = cat(1,XTest,TTest);
    UPred = forward(net,XTTest);

    % Calculate target.
    UTest = solveBurgers(extractdata(XTest),t,0.01/pi);

    % Calculate error.
    UPred = extractdata(UPred);
    err = norm(UPred - UTest) / norm(UTest);

    % Plot prediction.
    nexttile
    plot(XTest,UPred,"-",LineWidth=2);
    ylim([-1.1, 1.1])

    % Plot target.
    hold on
    plot(XTest, UTest,"--",LineWidth=2)
    hold off

    title("t = " + t + ", Error = " + gather(err));
end

legend(["Prediction" "Target"])

サポート関数

バーガース方程式の求解関数

関数 solveBurgers は、[3] で概説されているように、時間 t におけるバーガース方程式の真の解を返します。

function U = solveBurgers(X,t,nu)

% Define functions.
f = @(y) exp(-cos(pi*y)/(2*pi*nu));
g = @(y) exp(-(y.^2)/(4*nu*t));

% Initialize solutions.
U = zeros(size(X));

% Loop over x values.
for i = 1:numel(X)
    x = X(i);

    % Calculate the solutions using the integral function. The boundary
    % conditions in x = -1 and x = 1 are known, so leave 0 as they are
    % given by initialization of U.
    if abs(x) ~= 1
        fun = @(eta) sin(pi*(x-eta)) .* f(x-eta) .* g(eta);
        uxt = -integral(fun,-inf,inf);
        fun = @(eta) f(x-eta) .* g(eta);
        U(i) = uxt / integral(fun,-inf,inf);
    end
end

end

モデル損失関数

与えられた入力 $(x, t)$ について、ネットワークの出力 $u (x, t)$ がバーガース方程式、境界条件、および初期状態を満たすように強制することで、モデルに学習させます。特に、次の 2 つの量が損失の最小化に寄与します。

$loss = {MSE}_{f} + {MSE}_{u}$ ,

ここで、 ${MSE}_{f} = \frac{1}{N_{f}} \sum_{i = 1}^{N_{f}} {| f (x_{f}^{i}, t_{f}^{i}) |}^{2}$ および ${MSE}_{u} = \frac{1}{N_{u}} \sum_{i = 1}^{N_{u}} {| u (x_{u}^{i}, t_{u}^{i}) - u^{i} |}^{2}$ です。

ここで、 ${x_{u}^{i}, t_{u}^{i}}_{i = 1}^{N_{u}}$ は、計算領域の境界上の選点に対応し、境界条件と初期状態の両方を計上します。 ${x_{f}^{i}, t_{f}^{i}}_{i = 1}^{N_{f}}$ は、領域内部の点です。

${MSE}_{f}$ を計算するには、モデルの出力 $u$ の微分 $\frac{\partial u}{\partial t}, \frac{\partial u}{\partial x}, \frac{\partial^{2} u}{\partial x^{2}}$ を求める必要があります。

関数 modelLoss は、ネットワーク net、ネットワーク入力 X および T、初期状態および境界条件である X0、T0、および U0 を入力として受け取り、学習可能なパラメーターについての損失および損失の勾配を返します。

function [loss,gradients] = modelLoss(net,X,T,X0,T0,U0)

% Make predictions with the initial conditions.
XT = cat(1,X,T);
U = forward(net,XT);

% Calculate derivatives with respect to X and T.
gradientsU = dlgradient(sum(U,"all"),{X,T},EnableHigherDerivatives=true);
Ux = gradientsU{1};
Ut = gradientsU{2};

% Calculate second-order derivatives with respect to X.
Uxx = dlgradient(sum(Ux,"all"),X,EnableHigherDerivatives=true);

% Calculate mseF. Enforce Burger's equation.
f = Ut + U.*Ux - (0.01./pi).*Uxx;
zeroTarget = zeros(size(f),"like",f);
mseF = l2loss(f,zeroTarget);

% Calculate mseU. Enforce initial and boundary conditions.
XT0 = cat(1,X0,T0);
U0Pred = forward(net,XT0);
mseU = l2loss(U0Pred,U0);

% Calculated loss to be minimized by combining errors.
loss = mseF + mseU;

% Calculate gradients with respect to the learnable parameters.
gradients = dlgradient(loss,net.Learnables);

end

参考文献

Liu, Dong C., and Jorge Nocedal. "On the limited memory BFGS method for large scale optimization." "Mathematical programming" 45, no. 1 (1989): 503-528.
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis, Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations https://arxiv.org/abs/1711.10561
C. Basdevant, M. Deville, P. Haldenwang, J. Lacroix, J. Ouazzani, R. Peyret, P. Orlandi, A. Patera, Spectral and finite difference solutions of the Burgers equation, Computers & fluids 14 (1986) 23–41.