Using Neural Networks to Solve System of Ordinary Differential Equations(ODEs)
Advisor: Dr. Anuj Chaudhri
Internship Location: Google
Onsite Mentor: Bing Wang
Everything around us, from biology, stocks, physics, or even common life scenarios can be mathematically modeled by differential equations. From exponential growth/decay and simple harmonic motion equation to the Laplace equation, used to describe steady state heat conduction, and Navier-Stokes equation, used to describe fluid turbulent flow, differential equations can range from easy to solve to incredibly difficult. This project will be focused on solving systems of ordinary differential equations(ODEs). Traditional methods used to solve these systems of ODEs such as the Euler method have truncation and round-off errors as well as instability with stiff functions, functions with flat areas, that make it less accurate for higher order equations while methods such as Runge-Kutta, RDK4, is accurate to the fourth power but is 4 times as costly as the Euler method. This is where neural networks come in. Neural networks have an advantage over traditional methods because the non-linearity of these systems do not cause problems as the would for traditional methods, such as the stiff functions discussed above. Additionally, neural networks can handle high dimensionality with ease, which makes it perfect for large amounts of data with multiple interacting variables. Because neural networks are not slowed down by non-linearity and high dimensionality, they can predict systems of ODEs more accurately and faster than traditional methods. Through research and experimentation, I will create a neural network model using Python libraries to predict the Lorentz System, a chaotic system in which switching the initial conditions can drastically change the function. I chose this system because it is unstable and non-linear, which will be able to show the neural network’s capabilities to its fullest. I will compile the results of my experimentation, architecture of the final neural network, and comparison of efficiency to traditional methods into a research paper.
Project Blog
Week 0: Using Neural Networks to Solve Systems of Ordinary Differential Equations
Hello everybody, my name is Zheyu and welcome to my first blog!
My project is designing neural networks to solve systems of ordinary differential equations(ODEs). This project is aimed at using a different approach to systems of ODEs which may have an efficiency and accuracy factor times faster than that of traditional numerical algorithms, such as RDK4.
By testing multiple ways of modeling neural networks and comparing the resulting accuracy and efficiency, I strive to find which structures give the best solution. ODEs can be used to model many worldly phenomenon, and being able to solve them with greater accuracy and speed would lead to better predictions in events such as weather, fluid flow, population growth, spread of disease, and more. With many applications in the real world, neural networks can help supply society with more accurate predictions to act on. Better predictions of the spread of disease can lead to better policies. Better predictions of population growth can lead to more foresight in housing and infrastructure. More accurate modeling of fluids can lead to better fluid animation in video games and animation.
The first step towards my project is to learn and build my foundations in both linear algebra and neural networks. I have some books and courses lined up and I will start learning next week. I can’t wait to start the experimentation phase after I strengthen my foundations!
Week 1: Just the beginning...
Hello everyone!
Welcome back to my senior project blog. I spent this week mostly absorbing information and planning my initial approach to the problem. Let’s get straight into what I learned and why setup is crucial before jumping into experimentation!
Best Linear Algebra Online Course Ever
Important to neural networks is linear algebra, an area of math that I can’t seem to grasp no matter how many times I learn it. Linear algebra seemed like a bunch of random arbitrary rules that have no meaning. However, after going through Gilbert Strang’s MIT linear algebra lectures(highly recommend if you want to self learn), I found that there is a lot of intuition behind concepts like matrix multiplication, row echelon elimination, invertibility, linear combinations, spans, etc. He would refrain from stating rules until the intuition behind the rules are realized, which made learning a lot easier. On the website, there are even practice questions for students to do.
If anyone wants to follow these lectures, I recommend watching at 1.5x speed.
Reviewing Neural Networks
I went through the neural networks course I took on Coursera to remind myself of the structure of neural networks and the different kinds of activation functions such as ReLU, logistic regression, and softmax, and different loss functions like mean square error(not really used in practice), and binary cross entropy.
Setting up Github and virtual environments for jupyter notebook
Some set up before doing the project.
I will use Github to keep track of progress and virtual environments and jupyter notebook to set up the right python version and libraries for the project. This project will be coded using Pytorch.
Planned Approach
I will be solving the system of ODEs numerically. The neural network will be trained on a specific system, and will take in parameters as inputs and the missing variable(s) as the output.
Why this is important
The last time I went into a CS project without good planning, I had to reprogram everything on the last day. That was super painful…because I also had a presentation on my program that day. While simply guessing and checking all different combinations of neural network structures, activation functions, and loss functions might work, that won’t give me a good understanding on what the neural network is doing. If I don’t understand what I am doing, it will take me ages to figure anything out.
Thanks for reading and see you in the next blog!
Week 2: Getting the hang of Pytorch and Neural Networks
Hey folks, welcome back to my blog!
This week passed by so quickly, it felt like I wrote my last post only yesterday. Anyways, let’s get into what I have been working on this week!
The first few days, I spent learning linear algebra from both the Linear Algebra MIT Lectures as well as from the book Linear Algebra and Learning from Data. I especially focused on the properties of determinants, eigenvectors, and eigenvalues, which have been very fascinating! Determinants have 10 different properties, of which only 3 are needed to derive the rest. These first three properties are:
Determinant of the identity matrix is 1
If you exchange two rows of a matrix, you reverse the sign of its determinant from positive to negative or from negative to positive.
If we multiply one row of a matrix by t, the determinant is multiplied by t.
When deriving the rest, it gives a very intuitive understanding of what determinants are and how they utilize the symmetries of matrix operations. I won’t get into the rigorous details, but I am glad to finally understand the reasoning behind the seemingly random formula: ad - bc.
Eigenvalues and eigenvectors…I was first curious about them when my cosmology class touched upon it. It was hard to understand at first, but after following all the lectures, I have a better understanding of them. An eigenvalue and eigenvector is defined by: Ax = λx, where λ is the eigenvalue and x is the eigenvector. What this means is that the eigenvector, after being transformed by the matrix, will stay parallel to the original vector while being scaled by the eigenvalue. For a certain matrix, you can find its eigenvalues and eigenvectors by rewritten the above equation into: |Ax - Iλ| = 0, where I is the identity matrix. If the matrix is size n by n, the max number of eigenvectors and eigenvalues will be n. What I found the most interesting is that the sum of the eigenvalues will be equal to the diagonal of the matrix. A lot of times, eigenvectors and eigenvalues are actually really hard to find because eigenvalues can repeat…and there would be a missing eigenvector. Hopefully I won’t encounter them in my project…
Speaking of my project, let’s go over the direct progress I have made this week! Honestly, it’s not a lot because I ran into a lot of trouble while learning Pytorch and implementing neural networks, but that was expected. I went through the basic tensor concepts, nth degree matrixes, operations provided by the Pytorch library, such as flatten, matmul, stack, and more. After familiarizing myself with those, I read the implementation of a neural network. I liked it better than tensorflow because it separated activation functions and the neuron layers, which was easier to understand for me. To familiarize myself with these implementations, I tried to implement a convolutional neural network that can detect handwritten digits. It was a struggle, but in the end I got it to work. I definitely underestimated how long it takes to train a network… Next, I researched what neural network structures can work for ODEs, and decided on trying to use simple RNNs. I have a lot more research to do and a lot more progress to be made over the next few weeks!
That is all I have for you today! If you want to follow my progress, my Git Hub is linked here: https://github.com/Fox-King777/Senior-Project.git
Week 3: The struggle is real...
Hey folks, welcome back to my blog!
As you can read in the title, this week has been a struggle. After researching some neural networks to solve the preliminary problem of my project, modeling neural networks to solve simple ODEs, I found two great pytorch libraries.
First is Pytorch Lightning, which automates the training process in one class so that all I have to worry about is the actual model of the neural network. While this library makes it extremely easy to train neural networks, a lot of the process is thus abstracted away, which makes it harder to debug. While using this library, I have encountered tons of bugs and problems in terms of loading data and training data. I have yet to get used to the new library.
The next resource I found was Torchdyn, which is a github repository/library dedicated to Neural ODEs, which is a generalization to the Resnet architecture of neural networks. This type of neural network is depth invariant, and is structured in a way such that the depth of the neural network is continuous. Meaning, you can graph the progression of the neural network as it grows deeper. I researched two papers discussing what these neural ODEs are, and how they can be used: https://arxiv.org/abs/2008.02389, https://arxiv.org/abs/1806.07366. These research papers are actually quite new. At first I had a lot of trouble understanding the network structure and the paper in general, but eventually I started to understand the unique depth invariant nature of the structure. After gaining somewhat of an understanding, I followed a few of the tutorials the Torchdyn repository provided to learn how to implement Neural ODEs. Then I tried implementing this structure to model a simple exponential ODE, however the results were not ideal. Weirdly, as I graphed the results in terms of the depth of the neural network, I got an exponential looking graph. However, as I graphed the x to the predicted y, I got a function that is similar to the activation function I used. I have no idea how that occurred. It must be because I might have misunderstood or did not fully understand Neural ODEs. I tried again and again, tinkering with the structure, but nothing worked.
I took a step back, and tried using a MLP neural network with Pytorch Lightning, and still, the results did not match the function I was trying to model. This was very frustrating because this is only the preliminary problem, however, the good news is that because of the multidimensional nature of neural networks, generalizing this to a system of ODEs is a lot simpler than I thought it was. Once I get these neural networks to work with simple ODEs, I believe that this project will move forward a lot faster.
Although this week was less than ideal in terms of results, I feel that I am gaining a much better grasp of what I need to do to get this project finished. It is a lot more challenging than I anticipated, but also a lot more fascinating as well.
Thanks, and see ya guys next week!
Week 4: Abandoning NeuralODEs
Hey folks, welcome back to my humble blog…and this week…we have SUCCESS. I have a lot to share, so let’s get into it.
Because last week was such a struggle, I spent the entirety of Monday studying the basics of deep learning through the textbook Deep Learning by Ian Goodfellow, Yoshua Bengio, and Aaron Courville. This book is a fully comprehensive guide on deep learning, including linear algebra, probability, history and development of AI, information theory, and numerical optimization. Because of this, the textbook was super dense and it took a lot of time to read—but it definitely was definitely worthwhile as it helped me gain a better understanding of the neural network structure. Admitted, my other motivation for going through this textbook is because I wanted to learn more about deep learning.
Tuesday is the start of the trial and error process, and just like last week, it was tough. I tried to tinker with a pytorch lightning neural network model to fit an exponential equation but no matter what hyper parameters I used, the graph always looked like a flat line that suddenly asymptotes. Next, I decided to try to tinker with the NeuralODE model I implemented last week to fit a sine function but the graph the program produced was nothing close to a sine function. At this point, I was suspecting that NeuralODEs was not actually helping me accomplish my goals, so instead I used a lightning model to fit a sine function. This worked for short ranges only and did not have a high accuracy.
On Wednesday, I discovered that NeuralODEs were unbeneficial to my project goal. But I did not know that until my meeting with my external advisor. We figured out that NeuralODEs are simply a generalization of a residual network. The way residual networks solve ODEs is basically the same as using numerical methods to solve it. NeuralODEs are like residual networks except the depth of the neural network is modeled to be continuous instead of discrete. NeuralODEs are good at modeling data in which the correspondence between input and output is an unknown ODE, such as weather and fluids. This is not helpful to my project as I am trying to solve for a known and given system of ODEs. Though disappointing, it was able to point me in the right direction.
After running into many dead ends, I found this lecture: https://compphysics.github.io/MachineLearning/doc/pub/week43/html/week43-bs.html
With this resource, I was able to implement a neural network with one hidden layer to solve an exponential decay and growth differential. The ODE is given by dy/dx = x. This method differed from my past method because instead of generating data, it instead using a loss function that is modeled off of the given differential equation. The loss function is not calculated between the predicted y and the analytical y, but instead between the predicted y and the derivative, dy/dx. This derivative is calculated through calculating a gradient of the neural network. Because of the accuracy of this method, I plan on generalizing this method for a system of ODEs as well.
While there were many dead ends, I was finally able to implement a working program for a simple ODE. I don’t expect the rest of this project to be easy, however I have a clear plan as to how I am going to get it done.
That is all the progress I made this week, so see you guys next time!
Week 5: Time to move on
Hey folks, welcome back to my blog. This post will be a bit shorter than previous blogs because as of yesterday, I am in Taiwan. Nevertheless, I still made great progress so let’s get into it.
The program I implemented to solve exponential ODEs(using numpy autograd library) last week, which calculated loss between the estimated solution and gradient, was tested on trigonometry and polynomial functions. In both of these cases, the algorithm worked to a high degree of accuracy. For reproducibility, I tried to implement this model using Pytorch Lightning. And this is where I get stuck. I could not get this type of model to work in Pytorch because of the differences between Pytorch autograd and numpy autograd. It was difficult to debug because Pytorch Lightning automates the training process, abstracting parts of the code that might be causing a problem. The assumption that both autograd functions worked the same way was wishful thinking. I was comparing apples to oranges. Though I am not giving up on reimplementing it into Pytorch, I have to move on.
Using Pytorch Lightning, I implemented my initial approach of simulating x and y values and asking the model to fit the data. I failed to get this approach to work two weeks ago, however this week I was able to find the right hyperparameters for it to work. This model also worked with various functions to a high degree of accuracy, which I am happy about. The change I made to the hyperparameters was as simple as adjusting sample size, max epoch, and batch size. Additionally, the training process was also faster.
Moving forward with the project, I’ve decided to implement systems of ODEs in both of the ways I implemented above while reimplementing the autograd method in Pytorch Lightning on the side. Once I get a simple system of ODEs to work, I will move onto the Lorentz system.
That is all the progress I made this week, see you guys next time!
Week 6: Simple System of ODEs
Hey folks, welcome back to my blog. This week was a big mix of frustration and successes so let’s get into it.
On Monday, I attempted to debug my neural network model which calculated loss between the estimated solution and the gradient of the neural network, calculated using the pytorch auto gradient function. As I found out last week, this was difficult because my assumption that the pytorch autograd and numpy autograd worked the same way was inaccurate. I created a new jupyter notebook to try out all of the autograd features of both libraries on the same functions in order to find the difference. And it turns out, numpy autograd works on functions while pytorch autograd does not. Additionally, pytorch autograd works on pytorch functions such as the pytorch loss functions and the forward function for pytorch neural networks. However, if I multiply these functions by another variable, the entire thing does not work. For example, let’s say we have a forward function that takes in one parameter x, forward(x). Now let’s say we set a new variable A = x * forward(x). If we call autograd on A, this causes an error, which is a huge problem because that is the exact gradient I need to find to calculate the loss in each training step. I will need to find a work around.
On Tuesday, I attempted solving a simple system of ODEs: dy1/dx = y1, dy2/dx = y1 - y2, using the data fitting method. I created a pytorch data set which generated estimated y1 and y2 values iteratively using small steps in dx. After checking this generated data with the analytical solutions, I created a neural network with one input, x, and 2 outputs, y1 and y2. I had 5 hidden layers with the ELU activation function and used the MSE loss function. Then I started the training…and watched the loss after each epoch oscillate. That’s odd. This means that the loss was not converging. I stopped the training and lowered the learning rate to combat this but it did not work. I decided to just let it run for a few hundred epochs and check the results. Upon checking the results, I was surprised to find that the results are pretty accurate even though the MSE loss was in the thousands. It was close but not accurate enough for my liking. My guess is that there are a lot of local minimums in the loss function and the neural network is jumping between them, not able to find a suitable minimum. Tinkering with the depth and breadth of the neural network did not seem to help. I may need to try something else.
On Wednesday, I got sick of the slow training speed so I decided to set up some GPU acceleration. I downloaded CUDA and thought that I could automatically use GPU acceleration on my neural network. But NOPE. Tons of bugs. Man. First, pytorch was not even able to detect the gpu device. I tried to use the command prompt to troubleshoot but it turns out the command prompt couldn’t detect the python language. I found, painfully, that because I downloaded python from python.org instead of the Microsoft store, it was added to the path incorrectly. When I finally fixed that, I had to fix the issue where CUDA was not being detected by pytorch. I fixed this by redownloading pytorch with the CUDA computing option. GPU acceleration at first did not make the training process faster but once I tested it on a bigger neural network, I was definitely able to see a difference. I found that this is because GPU excels and parallel computing.
Because training was sped up, on Thursday I was able to tinker with the structure more efficiently. However I was not able to get any result from hours of trial and error. I think this is a dead end unless I find a significantly different structure. Frustrated, I decided to make my coding workspace more convenient by using VSCode instead of Jupyter Lab’s web compiler. I also decided to download Cygwin to learn Linux as that is good coding practice and a standard industry tool for file management and more.
Friday. HAPPY BIRTHDAY TO ME. Yep, today is my birthday. I am finally an adult. Maybe with my adult brain, I can be more successful with this project. Today I tried my data fitting with a second order system of ODEs, and realized that it doesn’t work. Which means either I am missing something, or that I need to switch to the autograd method. I think instead of getting stuck, I am going to try the autograd method.
That is all the progress I made this week, see you guys next time!
Week 7: Dead Ends and Black Magic
Hey folks, welcome back to my blog. As the title above suggests, I have hit an enormous dead end. This does not mean the project is over(I don’t give up that easily), but it does mean I have to try a different approach. Let’s go over how I hit the dead end and how I intend on getting out of it.
Monday was spent researching the properties of different activation functions, loss functions, and optimizers to find one that would fit my problem. I wanted to make the pytorch autograd model work. I altered the activation function of my model from sigmoid to RELU, leaky RELU, ELU, and GELU. While each activation was supposed to be an improvement from the previous one, it made almost no difference in the training process. Then I started tweaking the loss function, however I ended up back at MSE loss, my original loss function, as it was the most stable. Lastly, I tried changing the optimizers…and this is where I saw a big change. But unfortunately, this change was not good. Whether it was Adam, AdaMax, or Adadelta, they would always reach a local minimum loss that was super large. Although each behaved differently, none of them were able to converge near zero. Reaching a dead end here, I started a new notebook to try to figure out torch autograd, which I did. However, upon figuring it out, I am pretty sure my model should’ve worked.
On Tuesday, I tried altering the neural network layers to increase then decrease width taking inspiration from encoder and decoder structure. This did not change the training process at all. With nothing else to try, I gave up on torch autograd. Instead, I started training a neural network based on simulation to learn the Lorenz system. The result was not accurate even after running the training for 2 hours. On Wednesday, I tried tinkering with the model, but I was unable to increase the accuracy significantly. I think this is because the Lorenz system is simply just too difficult for the model to fit.
On Thursday, I felt incredibly demotivated. Pytorch just was not working for me. Both the simulation and autograd methods backfired. When I had my meeting with my external advisor, we went over what could have gone wrong but we were unable to gain a clear answer. After going over everything I have tried, we were genuinely stumped. Pytorch had hidden away too much of its complex functionality for me to effectively debug it. It felt like staring at...black magic. Strange and unexplainable. I am done dealing with Pytorch. I am going to code my own neural network without machine learning libraries.
Friday. I have a clear mathematical plan on how the algorithm works. By the end of the day, I implemented half of the neural network. I know exactly how I am going to finish this. But…I’ll leave my methods as a mystery for now. Y’all will have to tune in next week for when I reveal some magic of my own. MWAHAHAHA!
(This is how you do a cliffhanger right? Lol see you guys next week!)
Week 8: Almost There
Hey guys, welcome back to the blog. This week I have some exciting news so let's get into it.
As you know from last week, I decided to code the neural network from scratch instead of using machine learning libraries. A very daunting task at first but once I started coding, I started to realize I had an intuitive understanding of the neural network due to the first two weeks of preparation as well as the hours spent debugging pytorch code.
I approached the task by modularizing it, or in other words, splitting it into functions. I first coded the activation functions sigmoid, ELU, and tanh, then the MSE loss function. Afterwards I moved on the the harder tasks: initializing parameters, forward propagation, backward propagation, and gradient descent. Doing it this way lets me split the objective into smaller pieces which allows for much easier debugging. Debugging at each step is always the right way to go as it will avoid the trouble of running into black magic(reference from last week). I hate black magic. After running the training many times and figuring out how to match the dimensions of the matrices to perform forward and backward propagation, finally it worked. Well kind of. Randomly, I would get runtime overflow errors caused by NaN values, or Not a Number values. This means that somewhere in the math, the values go way too large. I also noticed that if I turn down the learning rate, the NaN values disappear. I reimplemented the parameters initialization and forward propagation more efficiently to see if that would make a difference. It did not fix the NaN values, but it improved the code and because of the way I organized the code, I didn’t have to change the other functions to fit the new changes.
At this point I was a little scared that this program was going to backfire but then I realized something as I was playing with the learning rates and observing the loss. As the learning rate increased, the loss was unable to converge. In fact, it diverged even faster. I realized that this is a common problem with neural networks and probably arises from the fact that I used a simple gradient descent algorithm. So I implemented the Adam optimizer which was meant to help with the rate of convergence and avoid saddle points in the loss function. And low and behold, it worked. My anxiety dropped at the same rate as the loss. Adam is a very interesting formula with physics inspired components. I encourage anyone who is interested in deep learning to read up on it. After graphing out the results, I still did not get the expected outcome but that is not the problem of the neural network but rather a problem with the trial solution. Trial solutions are assumptions on what the final function might look like. By adjusting it, I might just be able to find a solution to this problem.
Immense progress was made this week and I am learning a ton. Coding from scratch is a great way to learn things and I would suggest anyone who is interested in coding to do the same. I am excited to complete this project even if success is not guaranteed.
Thanks for reading my blog and see you next week.
Week 9: Black magic strikes once again
Hello everyone, welcome back to my blog. This week’s blog will be a bit shorter than usual because I have run into yet another roadblock with my code. Let’s get into the abnormalities I’ve observed and how I plan on dealing with them.
Picking up from last week, I noticed that although the loss converges, the neural network was not able to produce the right function. I concluded that the trial solution I based my neural network on was not the right form. So I began tinkering with the form of the equation, even adding trigonometric functions. This never worked, however, even with the loss converging. And this got me thinking. Theoretically, if the loss converges, that means the condition of the differential equation and the initial condition should be satisfied. If it could not be satisfied due to using an incorrect trial solution, the loss should not converge. This means that the initial trial solution was the right one. Then why did it not give the right solution to the system of differential equations?
I went back to the neural network I built to solve differential equations. I extended the range of the input. Big problem. Once I extended the input range from [0, 2] to [0, 10], it suddenly stopped working. Instead of giving me the expected graph, of exponential growth, it gave a graph similar to exponential decay. I tried many ranges, and [0, 2] was the only one that gave me the right solution. The others gave weird almost parabolic curves. This was perplexing and I still do not completely understand it. Another peculiar observation I made was that when I trained the model on [0, 2], it was somehow able to predict beyond the range it was trained on. Things are working and breaking in ways that they…should not be. I tried tempering with the ranges on the neural network model I used to solve the Lorentz system to see if I could get a similar result. I could not.
This could be a very huge issue and could prevent me from finishing this project if I am not able to figure out why these issues occur. However, I have plenty of things I could try. Either there is a logical problem in my code or there is a logical problem in the training of these models. I am first going to assume the latter and train in a mini-batch mode as opposed to the entire dataset at once. Then I am going to go through my code for logical bugs as well as print out the gradients to figure out what is happening. The model is falling into a trap somewhere and learning incorrectly. If somehow this doesn’t work, I’ll have to find another neural network model.
This weekend is going to be a lot of programming. What a difficult problem. Let’s do this.
Thanks for reading my blog and see you guys next week once I finish my final product!
Week 10: IT WORKS!
Hey guys, welcome back to my blog. Since I currently have a fever, I won’t be able to write much for this post but I’ll do my best to cover what I have accomplished.
Starting with fantastic news, my code works! Last week left me pretty demotivated and I really did not know if I was going to be able to get it working. After reading my code over and over again, I was unable to spot where it was going wrong. My code was so messy and unreadable that it was almost impossible to track whether it was math that was wrong or the algorithm that was wrong. I had to organize it. And that was what I spent my weekends doing. I modularized each part of the algorithm, testing each function to make sure none of the math is wrong. I wrote documentation for each of these functions, making sure I clearly understand the dimensions of the input and outputs. I made sure that these functions could work given a wide range of differential equations and trial solutions. After the implementation, I started testing.
Testing was a great success…until I tested the Lorenz system. The algorithm was unable to predict the Lorenz system. This is to be expected because the Lorenz system is a chaotic system. One small change in the initial state can create huge changes in a later state. This means if there is a small error in what the algorithm predicts, it can result in a graph that is entirely different than expected.
After testing the Lorenz system, I was going to implement the BFGS optimizer to see if I was able to get better results, however the day after I fell sick. Hopefully, I will recover enough to start testing some more and implement the BFGS optimizer during the weekends.
Thanks, this is all I have for this week's blog! See you guys at the Senior Project Presentations!