ОДУ первого порядка#
Общая информация#
В данном разделе мы на примере покажем, как применяется описанный метод в самом простом случае. Пусть стоит задача решения на отрезке \([a, b]\) уравнения
\[
\begin{equation*}
\frac{d f(x)}{d x}=G(x, f(x))
\end{equation*}
\]
с условием \(f(a) = A \in \mathbb{R} \).
Полагая (в вышеуказанных обозначениях) \(B \equiv A\), \(F(x, N(x ; \theta)) = (x-a) N(x ; \theta) \), получаем \(\hat{f} = A + (x-a) N(x ; \theta) \), то есть, что и требовалось, \(\hat{f}\) удовлетворяет граничному условию независимо от выхода нейросети.
Функция потерь для такого дифференциального уравнения будет задаваться выражением
\[
\begin{equation*}
J(\theta)=\frac{1}{M} \sum_{j=1}^{M}\left(\frac{\partial \hat{f}\left(x_{j} ; \theta\right)}{\partial x}-G\left(x_{j}, \hat{f}\left(x_{j} ; \theta\right)\right)\right)^{2}
\end{equation*}
\]
Также необходимо посчитать производную от пробной функции
\[
\begin{equation*}
\frac{\partial \hat{f}(x ; \theta)}{\partial x}=N(x ; \theta)+(x-a) \frac{\partial N(x ; \theta)}{\partial x}
\end{equation*}
\]
Пример 1#
Рассмотрим теперь конкретное уравнение (взятое из той же оригинальной статьи, указанной выше):
\[
\begin{equation*}
\frac{d f(x)}{d x}=x^{3}+2 x+x^{2} \frac{1+3 x^{2}}{1+x+x^{3}}-\left(x+\frac{1+3 x^{2}}{1+x+x^{3}}\right) f(x)
\end{equation*}
\]
с граничным условием \(f(0) = 1\)
Аналитическое решение (для сравнения результаты работы нейронной сети и истинного результата):
\[
\begin{equation*}
f(x)=\frac{e^{-\frac{x^{2}}{2}}}{1+x+x^{3}}+x^{2}
\end{equation*}
\]
Теперь мы готовы обучить простейшую нейронную сеть, используя Torch.
Наша нейронная сеть будет состоять из двух скрытых слоёв с регулируемым количеством нейронов. В качестве функции активации используем \(tanh\). Сначала попробуем обучить на 10000 эпохах с использованием Adam.
import torch
import torch.utils.data
import numpy as np
import matplotlib.pyplot as plt
import time
from torch.autograd import grad
class DataSet(torch.utils.data.Dataset):
"""
Класс для создания объектов, делающих torch-массивы в указанном промежутке
с указанным количеством точек
"""
def __init__(self, numSamples, xRange):
self.dataIn = torch.linspace(xRange[0], xRange[1], numSamples, requires_grad=True).view(-1,1)
def __len__(self):
return len(self.dataIn)
def __getitem__(self, idx):
return self.dataIn[idx]
class Fitter(torch.nn.Module):
"""
Собственно нейронная сеть, состоящая из одного входного и двух скрытых слоёв
с одинаковым количеством нейронов
"""
def __init__(self, numHiddenNodes):
super(Fitter, self).__init__()
self.fc1 = torch.nn.Linear(in_features = 1, out_features = numHiddenNodes)
self.fc2 = torch.nn.Linear(in_features = numHiddenNodes, out_features = numHiddenNodes)
self.fc3 = torch.nn.Linear(in_features = numHiddenNodes, out_features = 1)
def forward(self, x):
"""
Прямой ход нейронной сети
"""
h = torch.tanh(self.fc1(x))
h1 = torch.tanh(self.fc2(h))
y = self.fc3(h1)
return y
def plotNetwork(network, descentType, epoch):
'''
Функция для построения графиков и вывода результатов
'''
x = torch.linspace(xRange[0], xRange[1], 50, requires_grad=True).view(-1,1)
N = network.forward(x)
f_trial = trialFunc(x, N)
exact = solution(x)
MSECost = lossFn(f_trial, exact)
print("MSE между аналитическим и полученным решениями = ", MSECost.item())
exact = exact.detach().numpy()
x = x.detach().numpy()
N = N.detach().numpy()
plt.plot(x, trialFunc(x,N), 'r-', label = "Neural Network Output")
plt.plot(x, exact, 'b.', label = "True Solution")
plt.xlabel("x", fontsize = 16)
plt.ylabel("f(x)", fontsize = 16)
plt.legend(loc = "upper left", fontsize = 16)
plt.title(descentType + str(epoch) + " Epochs", fontsize = 16)
plt.show()
def trialFunc1(x, n_out):
"""
Тестовая функция. В вышенаписанных обозначениях B = 1, F(x, N) = x*N
"""
return 1 + x * n_out
def dTrialFunc1(x, n_out, dndx):
"""
Промзводная от тестовой функции
"""
return n_out + x * dndx
def diffEq1(x, f_trial, df_trial):
"""
По сути возвращает невязку при подстановке пробной функции"""
RHS = x**3 + 2*x + (x**2 * ((1+3*x**2) / (1 + x + x**3)))
LHS = df_trial + ((x + (1+3*(x**2)) / (1 + x + x**3) ) * f_trial)
return LHS - RHS
def solution1(x):
"""
Аналитическое решение задачи - для сравнения """
y = (torch.exp(-(x**2)/2) / (1 + x + x**3)) + x**2
return y
descentType = "Adam optimizer: "
def train(network, loader, lossFn, optimiser, numEpochs):
"""
Данная функция непосредственно обучает нейронную сеть. loader генерирует
батчи, lossFn - функция потерь
"""
cost_list=[]
network.train(True)
for epoch in range(numEpochs):
for batch in loader:
n_out = network.forward(batch)
dndx = grad(n_out, batch, torch.ones_like(n_out), retain_graph=True)[0]
f_trial = trialFunc(batch, n_out)
df_trial = dTrialFunc(batch, n_out, dndx)
diff_eq = diffEq(batch, f_trial, df_trial)
cost = lossFn(diff_eq, torch.zeros_like(diff_eq))
cost.backward()
optimiser.step()
optimiser.zero_grad()
cost_list.append(cost.detach().numpy())
network.train(False)
return cost_list
network = Fitter(numHiddenNodes=15)
xRange = [0, 1]
numSamples = 30
batchSize = 30
paramUpdates = 20000
train_set = DataSet(numSamples, xRange)
train_loader = torch.utils.data.DataLoader(dataset=train_set, batch_size=batchSize, shuffle=True)
lossFn = torch.nn.MSELoss()
optimiser = torch.optim.Adam(network.parameters(), lr=1e-3)
solution = solution1
trialFunc = trialFunc1
dTrialFunc = dTrialFunc1
diffEq = diffEq1
costList = []
totalEpochs = 10000
epoch = totalEpochs
start = time.time()
costList.extend(train(network, train_loader, lossFn, optimiser, totalEpochs))
end = time.time()
plotNetwork(network, descentType, epoch)
print(epoch, "final loss = ", costList[-1])
print("time = ", end - start, " s")
plt.semilogy(costList)
plt.xlabel("Epochs",fontsize = 16)
plt.ylabel("Cost",fontsize = 16)
plt.title("Training Cost",fontsize = 16)
plt.show()
MSE между аналитическим и полученным решениями = 4.019771626673219e-09
10000 final loss = 6.5465434e-07
time = 42.66474175453186 s
Попробуем для той же задачи другие встроенные инструменты: SGD и Rprop
descentType = "SGD optimizer: "
network = Fitter(numHiddenNodes=15)
xRange = [0, 1]
numSamples = 30
batchSize = 30
paramUpdates = 20000
train_set = DataSet(numSamples, xRange)
train_loader = torch.utils.data.DataLoader(dataset=train_set, batch_size=batchSize, shuffle=True)
lossFn = torch.nn.MSELoss()
optimiser = torch.optim.SGD(network.parameters(), lr=1e-3)
solution = solution1
trialFunc = trialFunc1
dTrialFunc = dTrialFunc1
diffEq = diffEq1
costList = []
totalEpochs = 10000
epoch = totalEpochs
start = time.time()
costList.extend(train(network, train_loader, lossFn, optimiser, totalEpochs))
end = time.time()
plotNetwork(network, descentType, epoch)
print(epoch, "final loss = ", costList[-1])
print("time = ", end - start, " s")
plt.semilogy(costList)
plt.xlabel("Epochs",fontsize = 16)
plt.ylabel("Cost",fontsize = 16)
plt.title("Training Cost",fontsize = 16)
plt.show()
MSE между аналитическим и полученным решениями = 7.396176079055294e-05
10000 final loss = 0.0031332911
time = 38.053231954574585 s
descentType = "Rprop optimizer: "
network = Fitter(numHiddenNodes=15)
xRange = [0, 1]
numSamples = 30
batchSize = 30
paramUpdates = 20000
train_set = DataSet(numSamples, xRange)
train_loader = torch.utils.data.DataLoader(dataset=train_set, batch_size=batchSize, shuffle=True)
lossFn = torch.nn.MSELoss()
optimiser = torch.optim.Rprop(network.parameters(), lr=1e-3)
solution = solution1
trialFunc = trialFunc1
dTrialFunc = dTrialFunc1
diffEq = diffEq1
costList = []
totalEpochs = 10000
epoch = totalEpochs
start = time.time()
costList.extend(train(network, train_loader, lossFn, optimiser, totalEpochs))
end = time.time()
plotNetwork(network, descentType, epoch)
print(epoch, "final loss = ", costList[-1])
print("time = ", end - start, " s")
plt.semilogy(costList)
plt.xlabel("Epochs",fontsize = 16)
plt.ylabel("Cost",fontsize = 16)
plt.title("Training Cost",fontsize = 16)
plt.show()
MSE между аналитическим и полученным решениями = 3.45977468896308e-08
10000 final loss = 9.24325e-06
time = 44.722307205200195 s
Эксперименты с размерами батчей и разными оптимайзерами показывают, что Adam является эффективным и универсальным решением, а в качестве размера батча лучше брать всю выборку (пока количество точек не слишком большое), выбирая достаточно большое количество эпох (не менее нескольких тысяч).
Пример 2#
Возьмём теперь другое уравнение первого порядка на отрезке \([0, 10]\)
\[
\begin{equation*}
\frac{d f(x)}{d x}=e^{-\frac{x}{5}} \cos (x)-\frac{1}{5} f(x)
\end{equation*}
\]
с условием \(f(0) = 0\).
Аналитическое решение
\[
\begin{equation*}
f(x)=e^{-\frac{x}{5}} \sin (x)
\end{equation*}
\]
Для его решения применим ровно тот же способ, что и ранее
def trialFunc1(x, n_out):
"""
Тестовая функция. В вышенаписанных обозначениях B = 0, F(x, N) = x*N
"""
return x * n_out
def dTrialFunc1(x, n_out, dndx):
"""
Промзводная от тестовой функции
"""
return n_out + x * dndx
def diffEq1(x, f_trial, df_trial):
"""
По сути возвращает невязку при подстановке пробной функции"""
RHS = df_trial + (1/5)*f_trial
LHS = torch.exp(-x/5) * torch.cos(x)
return LHS - RHS
def solution1(x):
"""
Аналитическое решение задачи - для сравнения """
y = torch.exp(-x/5) * torch.sin(x)
return y
descentType = "Adam optimizer: "
def train(network, loader, lossFn, optimiser, numEpochs):
"""
Данная функция непосредственно обучает нейронную сеть. loader генерирует
батчи, lossFn - функция потерь
"""
cost_list=[]
network.train(True)
for epoch in range(numEpochs):
for batch in loader:
n_out = network.forward(batch)
dndx = grad(n_out, batch, torch.ones_like(n_out), retain_graph=True)[0]
f_trial = trialFunc(batch, n_out)
df_trial = dTrialFunc(batch, n_out, dndx)
diff_eq = diffEq(batch, f_trial, df_trial)
cost = lossFn(diff_eq, torch.zeros_like(diff_eq))
cost.backward()
optimiser.step()
optimiser.zero_grad()
cost_list.append(cost.detach().numpy())
network.train(False)
return cost_list
network = Fitter(numHiddenNodes=15)
xRange = [0, 10]
numSamples = 50
batchSize = 50
paramUpdates = 20000
train_set = DataSet(numSamples, xRange)
train_loader = torch.utils.data.DataLoader(dataset=train_set, batch_size=batchSize, shuffle=True)
lossFn = torch.nn.MSELoss()
optimiser = torch.optim.Adam(network.parameters(), lr=1e-3)
solution = solution1
trialFunc = trialFunc1
dTrialFunc = dTrialFunc1
diffEq = diffEq1
costList = []
totalEpochs = 20000
epoch = totalEpochs
start = time.time()
costList.extend(train(network, train_loader, lossFn, optimiser, totalEpochs))
end = time.time()
plotNetwork(network, descentType, epoch)
print(epoch, "final loss = ", costList[-1])
print("time = ", end - start, " s")
plt.semilogy(costList)
plt.xlabel("Epochs",fontsize = 16)
plt.ylabel("Cost",fontsize = 16)
plt.title("Training Cost",fontsize = 16)
plt.show()
MSE между аналитическим и полученным решениями = 3.0525476176990196e-05
20000 final loss = 6.130902e-05
time = 94.68450498580933 s
