푸리에 신경 연산자를 사용하여 PDE 해 구하기

이 예제에서는 편미분 방정식(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는 다음과 같습니다.

$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \frac{1}{R e} \frac{\partial^{2} u}{\partial x^{2}}, x \in (0, 1), t \in (0, 1),$

여기서 $R e$ 는 레이놀즈 수입니다. 초기 조건은 다음과 같습니다.

$u (x, 0) = u_{0} (x), x \in (0, 1)$

이 방정식은 초기 조건 $u_{0}$ 을 사용하고 공간 영역 전체에 주기적 경계 조건을 적용하여 특정 시간에서의 해를 다음과 같이 정의합니다. $t = 1 : u_{0} ⟼ u (x, 1)$ .

초기 속도 $u_{0}$ 과 버거스 방정식의 해 $u (x, 1)$ 이 포함된 burgers_data_R10.mat에서 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

훈련 데이터 준비하기

이 예제에 지원 파일로 첨부된 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);

훈련 세트와 검증 세트의 크기를 확인합니다. 첫 번째 차원, 두 번째 차원, 세 번째 차원은 각각 배치 차원, 공간 차원, 채널 차원에 대응합니다.

size(U0Train)

ans = 1×3

        1638        1024           2

size(U0Validation)

ans = 1×3

         204        1024           2

신경망 아키텍처 정의하기

신경망 아키텍처를 정의합니다. 신경망은 여러 개의 푸리에 계층과 컨벌루션 계층으로 구성됩니다.

푸리에 계층 정의하기

주파수 영역에서 연산을 수행하는 신경망 계층을 만드는 함수를 정의합니다.

이 계층은 컨벌루션과 스펙트럼 컨벌루션을 적용하여 입력값을 처리한 다음 그 출력값을 합산합니다. 다음 도식은 이 계층의 아키텍처를 보여줍니다.

Diagram showing FNO layer architecture. The input is connected to both a spectral convolution layer and a convolution layer. The outputs of the two convolution layers are connected to an addition layer. The output of the FNO layer is the output of the addition layer.

스펙트럼 컨벌루션 계층은 주파수 영역에서 컨벌루션을 적용하며, 이를 통해 신경망이 복잡한 공간 종속성을 학습할 수 있기 때문에 PDE를 푸는 데 특히 유용합니다. 스펙트럼 컨벌루션 연산을 적용하기 위해 푸리에 계층은 이 예제에 지원 파일로 첨부되어 있는 사용자 지정 계층 spectralConvolution1dLayer를 사용합니다. 이 계층에 액세스하려면 예제를 라이브 스크립트로 여십시오.

function layer = fourierLayer(spatialWidth,numModes,args)

arguments
    spatialWidth
    numModes
    args.Name = ""
end
name = args.Name;

net = dlnetwork;

layers = [
    identityLayer(Name="in")
    spectralConvolution1dLayer(spatialWidth,numModes,Name="specConv")
    additionLayer(2,Name="add")];

net = addLayers(net,layers);

layer = convolution1dLayer(1,spatialWidth,Name="fc");
net = addLayers(net,layer);

net = connectLayers(net,"in","fc");
net = connectLayers(net,"fc","add/in2");

layer = networkLayer(net,Name=name);

end

신경망 만들기

신경망 아키텍처를 정의합니다. 이 예제의 신경망은 직렬로 연결된 여러 개의 Fourier-ReLU 블록으로 구성됩니다.

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.

다음과 같이 옵션을 지정합니다.

푸리에 계층에 대해, 공간 너비 16과 함께 모드 16개를 지정합니다.
크기가 [NaN spatialWidth 2]인 입력 계층을 사용하고 형식을 "BSC"(배치, 공간, 채널)로 지정합니다. 여기서 spatialWidth는 푸리에 계층의 공간 너비입니다.
convolution-ReLU 블록에 대해, 128개의 1×1 필터를 지정합니다.
마지막 컨벌루션 계층에 대해, 하나의 1×1 필터를 지정합니다.

numModes = 16;
spatialWidth = 64; 

layers = [
    inputLayer([NaN spatialWidth 2],"BSC");
    convolution1dLayer(1,spatialWidth,Name="fc0")
    fourierLayer(spatialWidth,numModes,Name="fourier1");
    reluLayer
    fourierLayer(spatialWidth,numModes,Name="fourier2")
    reluLayer
    fourierLayer(spatialWidth,numModes,Name="fourier3")
    reluLayer
    fourierLayer(spatialWidth,numModes,Name="fourier4")
    reluLayer
    convolution1dLayer(1,128)
    reluLayer
    convolution1dLayer(1,1)];

훈련 옵션

훈련 옵션을 지정합니다. 옵션 중에서 선택하려면 경험적 분석이 필요합니다. 실험을 실행하여 다양한 훈련 옵션 구성을 살펴보려면 실험 관리자 앱을 사용합니다.

SGDM 솔버를 사용하여 훈련시킵니다.
조각별 학습 계획을 사용하여 훈련시키고, 이때 초기 학습률을 0.001, 감소 인자를 0.5로 지정합니다.
크기가 20인 미니 배치로 Epoch 30회만큼 훈련시킵니다.
매 Epoch마다 데이터를 섞습니다.
입력 데이터의 차원이 배치 차원, 공간 차원, 채널 차원에 대응하므로 입력 데이터 형식을 "BSC"로 지정합니다.
훈련 진행 상황을 플롯에서 모니터링하고 검증 데이터를 지정합니다.
상세 출력값을 비활성화합니다.

schedule = piecewiseLearnRate(DropFactor=0.5);

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

상대 $L_{2}$ 손실 함수 정의하기

예측값과 목표값을 입력값으로 받고 상대 $L_{2}$ 손실을 반환하는 함수 relativeL2Loss를 만듭니다. 상대 $L_{2}$ 손실 함수는 $L_{2}$ 손실 함수의 변형으로, 예측값과 목표값 간 차이의 유클리드 노름을 목표값의 유클리드 노름으로 나누어 스케일을 반영합니다.

$R e l a t i v e L 2 (Y, T) = \sum_{n = 1}^{N} \frac{‖ T_{n} - Y_{n} ‖_{2}}{‖ T_{n} ‖_{2}}$

여기서

$Y$ 는 예측값을 나타냅니다.
$T$ 는 목표값을 나타냅니다.
$N$ 은 관측값 개수를 나타냅니다.

function loss = relativeL2Loss(Y,T)

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

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

loss = sum(p./q, 1);

end

신경망 훈련시키기

trainnet 함수를 사용하여 신경망을 훈련시킵니다. 사용자 지정 손실 함수 relativeL2Loss를 사용합니다. 기본적으로 trainnet 함수는 GPU를 사용할 수 있으면 GPU를 사용합니다. GPU를 사용하려면 Parallel Computing Toolbox™ 라이선스와 지원되는 GPU 장치가 필요합니다. 지원되는 장치에 대한 자세한 내용은 GPU 연산 요구 사항 (Parallel Computing Toolbox) 항목을 참조하십시오. GPU를 사용할 수 없는 경우, 함수는 CPU를 사용합니다. 실행 환경을 수동으로 선택하려면 ExecutionEnvironment 훈련 옵션을 사용하십시오.

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

신경망 테스트하기

홀드아웃 테스트 세트에 대한 예측값을 실제 레이블과 비교하고 평균제곱오차(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를 사용할 수 있으면 GPU를 사용합니다. GPU를 사용할 수 없는 경우, 함수는 CPU를 사용합니다. 실행 환경을 수동으로 선택하려면 testnet 함수의 ExecutionEnvironment 인수를 사용하십시오.

rmseTest = testnet(net,U0Test,U1Test,"rmse")

rmseTest = 
0.0112

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")

rmseTestZSSR = 
0.0113

새 데이터를 사용하여 예측하기

신경망을 사용하여 예측을 수행합니다.

새 데이터를 불러옵니다. 예시를 위해 처음 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를 사용할 수 있으면 GPU를 사용합니다. GPU를 사용할 수 없는 경우, 함수는 CPU를 사용합니다. 실행 환경을 수동으로 선택하려면 ExecutionEnvironment 옵션을 사용하십시오.

Y = minibatchpredict(net,U0New);

예측값을 플롯으로 시각화합니다.

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

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

참고 문헌

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.
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.

참고 항목

trainingOptions | trainnet | dlnetwork