{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "accelerator": "GPU", "colab": { "name": "AMGnotebook_ultimate_AF.ipynb", "provenance": [], "collapsed_sections": [] }, "kernelspec": { "display_name": "Python 3", "name": "python3" } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "uJeEfRup28rn" }, "source": [ "##An Investigation on the Optimal Parameters of the SOR and GMRES Solvers" ] }, { "cell_type": "markdown", "metadata": { "id": "ICgIn7JN8Sft" }, "source": [ "#0.1 Work Breakdown\n", "Class Taught by Professor Fournier and sponsored by Dr. Stephen Thomas\n", "\n", "Jacob Petersen: \n", "Data generation, building the AMG solver, computation of optimal strength of connection/relaxation parameter using fixed tolerance value and loops, optimal configuration of the post and pre-iterations for omega, and miscellaneous formatting and code optimization.\n", "\n", "Alex Semyonov: \n", "Computation of the theoretical optimal omega value, optimal configuration of the post and pre-iterations for theta, estimating theta using the standard form of a continuous optimization problem, miscellaneous formatting and code optimization.\n", "\n", "Christian Stout: \n", "Machine learning research and miscellaneous formatting and code optimization." ] }, { "cell_type": "markdown", "metadata": { "id": "SV2DyqcoB6Kv" }, "source": [ "#0.2 Abstract\n", "Algebraic multigrid (AMG) is an efficient method often utilized to solve large, sparse systems of equations. This method is advantageous relative to other multigrid methods as it does not require any geometric information and can be readily applied to various problems. The primary goal of this project is to determine the optimal parameters which enable the AMG to solve a system of equations with maximum efficiency. While working on this project, the team determined an optimal value for the relaxation parameter for the SOR algorithm and determined an optimal value for the strength of parameter for the GMRES solver." ] }, { "cell_type": "markdown", "metadata": { "id": "OdgoF4miB1C_" }, "source": [ "# 0.3 Introduction\n", "NREL (The National Renewable Energy Laboratory) is interested in the optimal construction of wind turbines. Flow dynamics are a particularly important factor to consider when constructing these wind turbines (as the air flow generated by a turbine may influence surrounding turbines). This air flow can be represented using the incompressible Navier Stokes equations which, when discretized, result in an exceptionally large system of equations. Solving these equations can often be very computationally taxing and time consuming. In order to efficiently solve this resulting system of equations, NREL has investigated the performance of several AMG methods in [“A Comparison of Classical and Aggregation-Based Algebraic Multigrid Preconditioners for High-Fidelity Simulation of Wind Turbine Incompressible Flows” by Thomas et al](https://www.osti.gov/pages/servlets/purl/1574702 ). This paper also examines the incompressible Navier-Stokes equations and their discretization. There are several parameters which may influence the computing speed and the generated residuals. By identifying the optimal parameters for a particular solving algorithm, one can reduce the number of iterations necessary to find a solution and minimize the residual associated with this solution.\n", "\n", "\n", " Ideally, the methods utilized to identify the optimal parameters can be generalized to other projects and perhaps can illuminate efficient methods to optimize the parameters of a given solver. As the discretization of partial differential equations often results in a large system of equations, identifying and optimizing the parameters which influence the iterations and residuals is an important task which may drastically influence computing requirements and solver accuracy. \n", " \n", " Our general approach to this optimization problem can be divided into four essential steps: \n", " \n", " 1. Generating the data (an *n x n* matrix **A** and an *n x 1* vector **b**). \n", " \n", " 2. Defining and building the appropriate solver (specifically one utilizing the SOR algorithm and another utilizing GMRES). \n", "\n", " 3. Optimizing the appropriate parameters using a brute-force approach (cycling through various parameter values and plotting residuals for each parameter)\n", "\n", " 4. Optimizing our method for determining the optimal parameter.\n", "\n", "When possible, the group computed optimal parameters both mathematically and numerically (although this was only possible with the SOR algorithm). This enabled the group to identify and explain potential discrepancies between the numerical and theoretical results. This, however, was more of a side-step than a progression step and thus is not included as an essential step. " ] }, { "cell_type": "markdown", "metadata": { "id": "aiaSLfPOB_zq" }, "source": [ "#1 Methods\n", "\n", "Generally, the method for each set of experiments use a `for` loop to vary the desired settings (omega, theta, pre-, and post-iterations) while constructing and running the AMG solver for each setting combination. \n", "\n", "`For` loops are used to construct and run the solver for each combination of settings--where the `for` loop is changing the desired setting, it mostly runs multiple values of omega or theta. If multiple nested for loops are used, it is to account for the combinations of pre- and post-smoother iteration settings and how those affect the solver as well. \n", "\n", "For the SOR solver, the theoretical optimal relaxation parameter was computed mathematically (using eigenvalues). For both the SOR and GMRES solvers, the optimal parameters ($\\theta$ and $\\omega$) are also calculated by minimizing an objective function with scipy.optimize.minimize_scalar. \n", "\n", "After computing an optimal parameter, various visual graphs are implemented to demonstrate how varying certain parameters influences the residuals and iteration counts.\n", "\n", "\n" ] }, { "cell_type": "code", "metadata": { "id": "Ba7PgFEDKuJq" }, "source": [ "# https://github.com/pyamg/pyamg\n", "!pip install pyamg" ], "execution_count": null, "outputs": [] }, { "cell_type": "code", "metadata": { "id": "AtdrJtWD3k82" }, "source": [ "import scipy\n", "import numpy\n", "import pyamg\n", "import pylab\n", "import random\n", "import pandas as pd \n", "import seaborn as sb\n", "import matplotlib.pyplot as plt\n", "from matplotlib.pyplot import figure\n", "from matplotlib import cm\n", "from random import seed\n", "from random import randint\n", "\n", "#Make sure this runtime is connected to a tf device\n", "#Edit -> Notebook Settings -> Hardware accelerator: GPU\n", "import tensorflow as tf\n", "from tensorflow import keras\n", "from keras import layers\n", "\n", "device_name = tf.test.gpu_device_name()\n", "if device_name != '/device:GPU:0':\n", " print('No tf GPU found')\n", " raise SystemError('GPU device not found')\n", "print(device_name)\n", "\n", "#Global Variables for consistancy in seed, size and tolerance level\n", "seednum= 5\n", "size = 100\n", "tol = 1e-13" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "Eumxk_n8GhTJ" }, "source": [ "##1.1 Problem Generation and AMG Construction\n", "\n", "After importing all of the necessary packages, the next step was to create matrices to run experiments on. To do this the `generate_problem` function was written, which produce a semi-random, sparse, and diagonally dominant matrix. This matrix is based off of a global seed variable to maintain consistency across all experiments. \n", "\n", "The next function is the `build_sor_amg` which is the function that is used when running experiments concering $\\omega$, where the $\\omega$ value and the pre- and post-smoother iterations can be set when the function is called. \n", "\n", "Finally, the last function, `build_theta_amg`, operates in the same way as previous function, except it focuses on the $\\theta$ parameter and allows the user to adjust it as necessary.\n" ] }, { "cell_type": "code", "metadata": { "id": "8NA-eyunVecd" }, "source": [ "def generate_problem(seednum, size, verbose=False):\n", " seed(seednum)\n", "\n", " diag = (randint(-size, size))+1000 # this is the random diagonal generated for the diagonal\n", "\n", " a = randint(-(abs(diag//6)),(abs(diag//6))) # these variables are the random stencil values\n", " b = randint(-(abs(diag//6)),(abs(diag//6)))\n", " c = randint(-(abs(diag//6)),(abs(diag//6)))\n", " if verbose:\n", " print(f'A= {a}\\nB= {b}\\nC= {c}\\nDiagonal= {diag}')\n", " \n", " stencil = [[a ,0 ,b], [c, diag, c], [b, 0, a]] # This stencil is set to be mostly sparse with its nonzero numbers generated randomly\n", " A = pyamg.gallery.stencil_grid(stencil, (size,size), dtype=float, format='csr')\n", " A = .5*(A + A.transpose()) # This makes A symmetric\n", " B = numpy.ones((A.shape[0],1))\n", " return A, B, diag\n", "\n", "def build_sor_amg(A, B, omega, preIts, postIts, maxCourse=10, verbose=False):\n", " ml = pyamg.smoothed_aggregation_solver(\n", " A, B, max_coarse=maxCourse, \n", " presmoother=('sor', {'sweep':'symmetric', 'omega':omega, 'iterations' : preIts}), \n", " postsmoother=('sor', {'sweep':'symmetric', 'omega':omega, 'iterations' : postIts}))\n", " \n", " if verbose:\n", " print (ml)\n", " b = numpy.random.rand(A.shape[0],1)\n", " x0 = numpy.random.rand(A.shape[0],1)\n", " return ml, b, x0\n", "\n", "def build_theta_amg(A, B, theta, preIts, postIts, maxCourse=10, verbose=False):\n", " ml = pyamg.smoothed_aggregation_solver(\n", " A, B, max_coarse=maxCourse,\n", " strength=('symmetric',{'theta':theta}), \n", " presmoother=('gmres', {'maxiter': preIts}), \n", " postsmoother=('gmres', {'maxiter' : postIts}))\n", " \n", " if verbose:\n", " print (ml)\n", " b = numpy.random.rand(A.shape[0],1)\n", " x0 = numpy.random.rand(A.shape[0],1)\n", " return ml, b, x0\n", "\n", "\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "S9CyB8908cW5" }, "source": [ "##1.2 Testing Omega with Default Pre and Post Iterations Values\n", "\n", "Using the two functions that are defined above, the problem matrix is generated and the SOR-based AMG is constructed. Then using a command from the PyAMG package we can just run the solver on our defined problem. It is in this step that we are varying $\\omega$, the relaxation parameter for SOR. The optimization of $\\omega$ is the first step in optimizing the efficiency of the AMG solver with SOR pre- and post-smoothing methods. \n", "\n", "\n", "The program starts $\\omega$ at 0.05 and ends at 2.0, testing every value of $\\omega$ with a step size of 0.025. \n", "In these tests, the values for the pre- and post-iterations are set to their default values of 1. This is just to analyze how much a change in $\\omega$ can influence the resulting residuals and iteration counts. \n", "\n", "\n" ] }, { "cell_type": "code", "metadata": { "id": "HGS3jxblbyHp" }, "source": [ "#This cell may take some time (typically under 30 seconds) to run\n", "omega_experiments = []\n", "A, B, diag= generate_problem(seednum, size)\n", "for omega in numpy.linspace(0.05, 2.0, 80):\n", " ml, b, x0 = build_sor_amg(A, B, omega, 1, 1)\n", " \n", " residuals = []\n", " x = ml.solve(b=b, x0=x0, tol=tol, residuals=residuals) \n", "\n", " omega_experiments.append((omega, residuals))\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "wkUU9phKfLdh" }, "source": [ "#2 Computing the Optimal Parameters" ] }, { "cell_type": "markdown", "metadata": { "id": "0cYrZLzDL4RW" }, "source": [ "## 2.1 Calculating the Theoretical Optimal $\\omega$\n", "When working with the SOR algorithm, $\\omega$ denotes the relaxation parameter. This relaxation parameter can drastically influence convergence rate. $\\omega$ can be mathematically calculated by finding the [spectral radius](https://en.wikipedia.org/wiki/Spectral_radius#Matrices) of the Jacobi iteration matrix and utilizing a [formula](https://en.wikipedia.org/wiki/Successive_over-relaxation#Convergence_Rate ) which allowed us to directly compute the optimal relaxation parameter. There are several assumptions made in this calculation, namely: (1) The Jacobi Iteration Matrix has real eigenvalues, (2) that Jacobi's method is convergent, and that a unique solution exist. Furthermore, our eigenvalue computation relies on A being symmetric. If A is not symmetric, the eigenvalue computation can be changed to eigs() but will take longer to compute." ] }, { "cell_type": "code", "metadata": { "id": "5IJQZtEIUcXM" }, "source": [ "\n", "A, B, diag = generate_problem(seednum, size)\n", "jacobi=scipy.sparse.identity(size*size)-A/diag #Definition of Jacobi matrix\n", "max_eigenval =scipy.sparse.linalg.eigsh(jacobi, k=1, which='LM',return_eigenvectors=False)#Returns Eigenvalue of the largest magnitude\n", "print(\"Largest Magnitude Eigenvalue: \", max_eigenval)\n", "optimal_omega = 1 + (max_eigenval/(1+numpy.sqrt(1-max_eigenval**2)))**2#Formula for optimal relaxation parameter\n", "print(\"Optimal Omega: \", optimal_omega)\n", "\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "5-nHK1KSiboL" }, "source": [ "##2.2 Continuous Optimization of $\\omega$\n", "\n", "The following code block finds the optimal $\\omega$ by first defining an [objective function](https://en.wikipedia.org/wiki/Optimization_problem#Continuous_optimization_problem): $\\omega$2iter($\\omega$) which returns the iteration count at which the residual falls below the desired tolerance level. The function then finds the $\\omega$ that minimizes the iteration count. Comparing the produced value to the theoretical omega demonstrates that the difference between the two values is rather miniscule, which suggests that the algorithm can accurately estimate the optimal relaxation parameter. " ] }, { "cell_type": "code", "metadata": { "id": "sYhZ6hp8iQpu" }, "source": [ "#This code was made by prof Fournier, I posted it here to try to disect it and compare with theoretical omega below.\n", "from scipy.optimize import minimize_scalar\n", "A, B, diag= generate_problem(seednum, size)\n", "#This is our objective function\n", "def ω2iter(ω) :\n", " ml, b, x0 = build_sor_amg(A, B, ω, 1, 1)\n", " residuals = []\n", " x = ml.solve(b=b, x0=x0, tol=tol, residuals=residuals)\n", " iterCount = numpy.nonzero([r == 0#< tolerance is set to 0 \n", " for r in residuals])[0] # 0th of all indexes of all residuals == 0\n", " if iterCount.size > 0 : \n", " return float(iterCount) \n", " else :\n", " return float(len(residuals))\n", "res_omega = minimize_scalar(ω2iter, bounds=(0.05, 2.0), method='Bounded', tol=0) #finds ω that minimizes number of iterations to cross desired tolerance\n", "res_omega" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "6z4k-5w0RVdr" }, "source": [ "##2.3 Using the Theoretical Optimal $\\omega$ To Find the Best Iteration Settings\n", "\n", "In this code, we have fixed the relaxation parameter to be the theoretical optimal parameter calculated above. This allows us to identify which pre/post-iteration settings are optimal in terms of speed and accuracy. To do this, the AMG construction and solver are repeated 25 times for all possible combinations of pre- and post-iterations while keeping the $\\omega$ value constant between experiments. \n" ] }, { "cell_type": "code", "metadata": { "id": "lXGI4He0MIIe" }, "source": [ "iteration_experiments = []\n", "A, B, diag= generate_problem(seednum, size)\n", "for preIts in numpy.arange(1, 6, 1): #Pre iterations loop\n", " for postIts in numpy.arange(1, 6, 1): #Post iterations loop\n", " ml, b, x0 = build_sor_amg(A, B, optimal_omega, preIts, postIts) #Optimal omega from above cell\n", "\n", " residuals = []\n", " x = ml.solve(b=b, x0=x0, tol=tol, residuals=residuals) \n", "\n", " iteration_experiments.append((preIts, postIts, residuals))" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "FKwSx7TERA2F" }, "source": [ "##2.4 Testing to find the Optimal $\\theta$\n", "\n", "The optimal $\\theta$ for running the solver cannot be calculated theoretically in a similar manner to the optimal $\\omega$.\n", "\n", "Instead the use of a brute force method is applied, in which a number of experiments are done for each combination of pre- and post-smoother settings. There are 25 different combinations for the pre- and post-smoother settings, which can have an integer value from 1 to 5 each, and the range of $\\theta$ spans 0.0 to 0.6. For every comination of pre- and post- iterations, the program will construct and solve an AMG for every $\\theta$ value from 0.0 to 0.6 with steps of 0.01. \n", "\n", "After each solution is finished, the residuals list is passed to a new array to store the values along with the current $\\theta$ and pre- and post-iterations. \n", "\n" ] }, { "cell_type": "code", "metadata": { "id": "gAh0K5Eyi3AI" }, "source": [ "# THIS TESTS OVER 1500 SETTINGS COMBINATIONS, TAKES A WHILE TO RUN#\n", "\n", "theta_experiments = []\n", "hm_data = []\n", "hm_data_all = []\n", "seednum= 5\n", "A, B, diag = generate_problem(seednum, size)\n", "\n", "for preIts in numpy.arange(1, 6, 1): # Pre iteration loop\n", " for postIts in numpy.arange(1, 6, 1): # Post iteration loop\n", "\n", " temp_resid_len = []\n", " temp_theta_store = []\n", "\n", " for theta in numpy.linspace(0.01, .60, 60): # Theta value loop\n", " ml, b, x0 = build_theta_amg(A, B, theta, preIts, postIts)\n", "\n", " residuals = []\n", " x = ml.solve(b=b, x0=x0, tol=tol, residuals=residuals) \n", "\n", " # Array stores ALL data for display later\n", " theta_experiments.append((theta, preIts, postIts, residuals))\n", " \n", " # Temporary storage for each combination of pre and post its\n", " temp_theta_store.append((theta, preIts, postIts, residuals/residuals[0]))\n", "\n", " # Temporary storage of length for all thetas in current combination of pre and post its \n", " temp_resid_len.append(len(residuals)) \n", " \n", " # These two lines find the theta with the fewest iterations for each pre and post it combination\n", " minpos = temp_resid_len.index(min(temp_resid_len))\n", " the, pre, post, resid = temp_theta_store[minpos]\n", "\n", " # Storage for the 25 optimal thetas and their respective results for use in displaying as heatmaps\n", " hm_data_all.append((the, pre, post, resid[-1]/resid[0], len(resid)))\n", " hm_data.append((the, pre, post))" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "GjCooyu0isrl" }, "source": [ "##2.5 Continuous Optimization Problem $\\theta$\n", "This code block once again defines an objective function: Θiter(Θ). Which is a function of $\\theta$ that returns the iteration count (to hit desired tolerance). This function is then minimized using scipy.optimize.minimize_scalar(). The produced theta solution always tends towards the upper bound and differs from the estimated theta value discovered by utilizing the \"brute-force\" approach above." ] }, { "cell_type": "code", "metadata": { "id": "rz22EUSfiew2" }, "source": [ "\n", "from scipy.optimize import minimize_scalar\n", "A, B, diag = generate_problem(seednum,size)\n", "def Θiter(Θ): #Objective function definition\n", " \n", " m1, b, x0 = build_theta_amg(A,B,Θ,1,1) #Varying pre/post iteration essentially has no effect on the solution\n", " residuals = []\n", " x = ml.solve(b=b,x0=x0,tol=tol,residuals=residuals)\n", " iter_count = numpy.nonzero([r < tol for r in residuals])[0]\n", " if iter_count.size > 0 :\n", " return float(iter_count) \n", " else :\n", " return float(len(residuals))\n", "\n", "res_omega = minimize_scalar(Θiter, bounds = (0.01, 0.6), method = 'Bounded', tol=0) #The smaller the tolerance the closer res_omega is to upper bound. Theta is between 0 and 0.60\n", "res_omega" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "C2iP3z3siBv_" }, "source": [ "##2.6 Pre- and Post- iterations with a Fixed $\\theta$\n", "\n", "The section above tests all values for theta against every combination of pre- and post-iterations, but that also takes a long time to run. This section allows the user to run an individually chosen theta value on all combinations of pre- and post-iteration values. The main reason we left this in is to test any value that might be of interest, or to test fewer values by hand, which will result in much faster response times for ease of use. \n" ] }, { "cell_type": "code", "metadata": { "id": "nvW-EJDHj6nD" }, "source": [ "theta_iteration_experiments = []\n", "\n", "\n", "theta=0.10\n", "A, B, diag= generate_problem(seednum, size)\n", "for preIts in numpy.arange(1, 6, 1):\n", " for postIts in numpy.arange(1, 6, 1):\n", " ml, b, x0 = build_theta_amg(A, B, theta, preIts, postIts)\n", " residuals = []\n", " x = ml.solve(b=b, x0=x0, tol=tol, residuals=residuals) \n", "\n", " theta_iteration_experiments.append((preIts, postIts, residuals))" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "fLCxmHmeeRdS" }, "source": [ "#3 Results\n", "\n", "Seperated into different categories below are the results, and their graphs, for each experiment run above. " ] }, { "cell_type": "markdown", "metadata": { "id": "g0UbdJWN81tG" }, "source": [ "##3.1 Testing Omega with Default Pre and Post iterations Values: Results\n", "Once the AMG solver has run for all of the desired values of $\\omega$, all of the resulting residuals will be graphed according to their total number of iterations for a solution. \n", "\n", "Below the graph, the optimal omega value from all of the tested values is printed. In this case, the optimal $\\omega$ was determined by sorting the list of residual arrays by the length of the residual arrays. This gives the $\\omega$ value that corresponds to the fewest number of residuals, or fewest iterations. \n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "code", "metadata": { "id": "TJ7_ku55MZfp" }, "source": [ "residuals_len =[]\n", "figure(figsize=(18, 12))\n", "plt.xlim(0,20)\n", "\n", "# Loop to graph each residuals list from the experiments\n", "for experiment in omega_experiments:\n", " omega, residuals = experiment\n", " #Fournier added the edits below, the color is sorted, fewest iterations are light green and become more purple/red with more iterations\n", " pylab.semilogy(residuals/residuals[0], '.-', label=f\"ω: {omega:.2f}\",\n", " c=cm.rainbow((omega - omega_experiments[0][0])/(omega_experiments[-1][0] - omega_experiments[0][0]))) \n", " residuals_len.append(len(residuals/residuals[0]))\n", "\n", "# Graph Labels\n", "pylab.xlabel('iterations')\n", "pylab.ylabel('normalized residual')\n", "pylab.legend(loc='upper right',ncol=4)\n", "pylab.show()\n", "\n", "# Just some useful information as well as the optimal choice for Omega in regards to iterations\n", "minpos = residuals_len.index(min(residuals_len))\n", "print(f'Lowest iteration count = {min(residuals_len)} \\n')\n", "print(f'All iteration counts = {residuals_len} \\n')\n", "print(f'Optimal Omega index = {minpos} \\n')\n", "print(f'Optimal Omega value = {omega_experiments[minpos]}')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "2vfI6bQ8HW3h" }, "source": [ "##3.2 Using Theoretical Optimal $\\omega$ to find Iteration Settings: Results\n", "\n", "The graph below suggests that there is a connection between the pre- and post-smoother iterations, and how many overall iterations it takes for the AMG solver to finish. It seems that as the pre- and post-iterations both increase, the number of overall iterations decreases. Unfortunately, a problem was identified in the code block below just prior to the end of the term and the group was unable to correct it in time. Changing the direction of the for-loop, that is, looping through iteration_experiments in the reverse order produces a different result. This is troubling as the answer produced should be the same regardless of the order of the loop and must be investigated further." ] }, { "cell_type": "code", "metadata": { "id": "MgtJ4DqaMziE" }, "source": [ "residuals_len =[]\n", "figure(figsize=(20, 16))\n", "plt.xlim(0,20)\n", "\n", "for experiment in iteration_experiments:\n", " preIts, postIts, residuals = experiment\n", " residuals_len.append(len(residuals))\n", " pylab.semilogy(residuals/residuals[0], '.-', label=f\"Pre Its: {preIts} Post Its: {postIts}\",\n", " c=cm.rainbow((preIts-iteration_experiments[0][0])/(iteration_experiments[-1][0]-iteration_experiments[0][0])), #once again the prof switched up the color scheme to make the graph more presentable\n", " marker='${}$'.format(postIts))\n", " \n", "pylab.xlabel('iterations')\n", "pylab.ylabel('normalized residual')\n", "pylab.legend(loc='best',ncol=2)\n", "pylab.show()\n", "\n", "minpos = residuals_len.index(min(residuals_len))\n", "print(f'Lowest iteration count = {min(residuals_len)} \\n')\n", "print(f'All iteration counts = {residuals_len} \\n')\n", "print(f'Optimal Iteration Index = {minpos} \\n')\n", "print(f'Optimal Iteration Values = {iteration_experiments[minpos]}')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "tnyTIHf3RWbB" }, "source": [ "##3.3 All $\\theta$ Values Tested for all Pre and Post Iteration Combinations: Results\n", "\n", "The graph below shows the results from the experiments in the prior cell. It is graphed using the same method described previously, with one caveat. Despite best efforts, removing the repeating final residual from many of these experiments would also remove many data points that are important. Any method attempted to remove these values also removes residual values that are below the desired tolerance level, sometimes leaving only values that are above the tolerance level. \n", "\n", "The values of $\\theta$ and pre and post-iterations that led to the repeating of iterations, despite having already hit the tolerance threshold, are quite troublesome. It seems that some part of the AMG package does not terminate the algorithm sometimes, even though it seems to hit the desired tolerance level prior to repeating the last residual a hundred times before stopping. This leads to another problem that came up, there is no actual guarantee in these experiments that the final residual achieved is below the tolerance level. The PyAMG package is designed to solve systems with this tolerance constraint, but there was never a check implemented to make sure this is the case. That is an oversight of the project and something to keep in mind when viewing the below results.\n", "\n", "While the method for finding the best $\\theta$ using the residuals list with the shortest length still functions, some part of the PyAMG package causes the solver to repeat iterations of the last residual. So for the sake of accuracy with the existing values, the residuals list is kept the same. " ] }, { "cell_type": "code", "metadata": { "id": "mbXDZ8ZIi30i" }, "source": [ "residuals_len = []\n", "figure(figsize=(20, 16))\n", "plt.xlim(0,20)\n", "\n", "# Loop to graph each residuals list from the experiments\n", "for experiment in theta_experiments:\n", " theta, preIts, postIts, residuals = experiment\n", "\n", " # Commented out the legend because it was too long\n", " pylab.semilogy(\n", " residuals, '.-', label=f\"θ: {theta:.2f}\",c=cm.rainbow((theta-theta_experiments[0][0])/ # Fournier\n", " (theta_experiments[-1][0]-theta_experiments[0][0])))\n", " \n", " residuals_len.append(len(residuals))\n", "\n", "pylab.xlabel('Iterations')\n", "pylab.ylabel('Normalized Residual')\n", "#pylab.legend(loc='best')\n", "pylab.show()\n", "\n", "# Prints the Theta and Pre/Post iterations combo that used the least iterations to solve\n", "minpos = residuals_len.index(min(residuals_len))\n", "print(f'Optimal Theta index = {minpos} \\n')\n", "print(f'Optimal Theta value = {theta_experiments[minpos]}')\n", "print(min(residuals_len))" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "Yb2p9HHTzE5Y" }, "source": [ "##3.4 Heatmap View of the Optimal $\\theta$ for Each Combination of Pre- and Post-Iterations\n", "\n", "First, a small note about the heatmaps, the axes are indexed from 0 to 4 even though the pre and post iteration values are from 1 to 5, so keep that in mind when viewing the results. For example, if you want 3 pre-iterations and 4 post-iterations, look at row 2 column 3.\n", "\n", "The heatmap on the left shows the $\\theta$ value that resulted in the fewest number of iterations for each combination of pre- and post-iterations. Most of these were 0.01, but as the map shows there were a few values that differed.\n", "\n", "The second heatmap shows the final residual when the solver is run with the $\\theta$ settings from the first heatmap. \n", "\n", "The final heatmap shows how many iterations it took to reach the final residuals shown in the second heatmap. The problem with these heatmaps lies in the third and final one because of the combinations that led to over 100 iterations. These are not indicative of a proper solution by the solver and so they skew the results shown in the heatmap, since they are not a reliable indicator of how that combination of $\\theta$ and pre- and post-iterations truly affect the solution. There is a possibility that these solutions might be important but the package seems to repeat the residual without finishing the solutions, and we could not find out why. \n", "\n" ] }, { "cell_type": "code", "metadata": { "id": "-n6RNFshy_DY" }, "source": [ "theta_map = numpy.zeros((5,5))\n", "resid_map = numpy.zeros((5,5))\n", "its_map = numpy.zeros((5,5))\n", "k = 0\n", "\n", "#for x in hm_data_all:\n", " #print(x)\n", "\n", "for i in numpy.arange(0,5,1):\n", " for j in numpy.arange(0,5,1):\n", " the, pre, post, resid, len_resid = hm_data_all[k]\n", "\n", " #print(i, j)\n", " \n", " theta_map[i][j]= the\n", " resid_map[i][j]= resid\n", " its_map[i][j]= len_resid\n", "\n", " k += 1\n", "\n", "fig, ax = plt.subplots(ncols=3, sharey=False)\n", "fig.set_size_inches(24,6)\n", "\n", "\n", "the_heat = sb.heatmap(theta_map, annot= True, ax=ax[0], linewidths= 0.5, cmap=\"YlGnBu\")\n", "resid_heat = sb.heatmap(resid_map, annot=True, ax=ax[1], linewidths= 0.5, cmap=\"YlGnBu\")\n", "its_heat = sb.heatmap(its_map, annot=True, ax=ax[2], linewidths= 0.5, cmap=\"YlGnBu\")\n", "\n", "fig.suptitle('Optimal Theta Results for All Combinations of Pre- and Post-Iterations')\n", "ax[0].set_title('Optimal Theta')\n", "ax[1].set_title('Last Residual')\n", "ax[2].set_title('Total Iterations')\n", "\n", "for ax in ax.flat:\n", " ax.set(xlabel='Post-Iterations', ylabel='Pre-Iterations')" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "dcHjY9l0WTZG" }, "source": [ "##3.5 Pre and Post iterations with a Fixed $\\theta$: Results\n", "\n", "This code block below graphs the results of fixing a $\\theta$ value and varying the pre- and post-iterations, in the same manner as the previous graphing cells. Excluding the results cell for section *Testing to Find the Optimal $\\theta$*. \n", "\n", "The results once again show that there is a relationship between the solver iterations and the number of iterations the pre- and post-smoothers have to run. If the pre- and post-smoothers are allowed more iterations each, then the overall iteration count of the solver goes down. \n", "\n", "This is the case if the repeating residuals are excluded, as mentioned above, the solver sometimes continues running even after the desired tolerance is reached. The group was unable to identify the cause for this troubling finding.\n" ] }, { "cell_type": "code", "metadata": { "id": "8euM5ZMRUy69" }, "source": [ "\n", "residuals_len =[]\n", "figure(figsize=(20, 16))\n", "plt.xlim(0,20)\n", "\n", "\n", "for experiment in theta_iteration_experiments:\n", " preIts, postIts, residuals = experiment\n", " residuals_len.append(len(residuals))\n", " pylab.semilogy(residuals/residuals[0], '.-', label=f\"Pre Its: {preIts} Post Its: {postIts}\",\n", " c=cm.rainbow((preIts-theta_iteration_experiments[0][0])/(theta_iteration_experiments[-1][0]-theta_iteration_experiments[0][0])),\n", " marker='${}$'.format(postIts))\n", " \n", "\n", "pylab.xlabel('iterations')\n", "pylab.ylabel('normalized residual')\n", "pylab.legend(loc='best', ncol=2)\n", "pylab.show()\n", "\n", "\n", "minpos = residuals_len.index(min(residuals_len))\n", "print(f'Lowest iteration count = {min(residuals_len)} \\n')\n", "print(f'All iteration counts = {residuals_len} \\n')\n", "print(f'Optimal Iteration Index = {minpos} \\n')\n", "print(f'Optimal Iteration Values = {iteration_experiments[minpos]}')\n", "print(theta)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "dP-HIL73yjrj" }, "source": [ "##3.6 Discussion \n", "\n", "The difference between the tested optimal $\\omega$ and the theoretical optimal $\\omega$ is the most interesting result. This result was unexpected, and raised a lot questions as to why the tested $\\omega$ and the theoretical $\\omega$ were always different. Looking a bit more into the mathematics of iterative methods, it seems that the optimal theoretical $\\omega$ is considered optimal in terms of the rate of convergence of the iterative method. The experiments regarding $\\omega$ revealed the value that reduces the number of iterations for the solver. This does not imply that this value also acts as the optimal $\\omega$ for the rate of convergence.\n", "\n", "The tests involving the strength of connection paramter, $\\theta$, also had some interesting results, each experiment showed that there was not a heavy dependence on $\\theta$ for the solver to run to completion. Any value used seemed to have the same effect on the sparse, diagonally dominant matrices. This might not be the case if given matrices with different properties, which is something to look into for further testing. Generally, increasing both the pre- and post- iterations decreases the solver iterations; this result was achieved regardless of the values of the other parameters ($\\omega$, $\\theta$).\n", "\n", "The strength of connection parameter ($\\theta$) was also optimized by minimizing a function that takes in $\\theta$ and returns the iteration counts required to reach a desired tolerance level. While this approach worked wonderfully for optimizing the relaxation parameter, $\\omega$, used in the SOR solver, using it to compute an optimal theta yielded unclear results. The supposed optimal theta value always tends to the upper bound. So for example if the upper bound is set to 0.6, the optimal theta will supposedly be 0.5999999, if the upper bound is set to 0.3 it will 0.299999. This occured regardless of what post/pre-iteration settings were used. While the reason for this was never unequivocally determined it may be due to the construction of the GMRES solver in the PyAMG source code and should be investigated in the future. Since this process accurately computed the optimal relaxation parameter, the group is inclined to believe that the issues lies with the GMRES solver and not the minimize_scalar() function." ] }, { "cell_type": "markdown", "metadata": { "id": "T_mVdz60uv5O" }, "source": [ "#4 Conclusions" ] }, { "cell_type": "markdown", "metadata": { "id": "k4or6ec5u_-e" }, "source": [ "\n", "\n", " Due to restrictions present in Google Colaboratory, the group was not able to work directly with \n", "the data produced by NREL. This was simply because the data produced by NREL was too large for Colab to handle. After collaborating with both Dr. Aimé Fournier and Dr. Stephen Thomas, the group successfully generated data (a diagonally dominant A matrix and **b** vector for the solver to ingest) and identified which parameters appeared to influence the number of iterations as well as the generated residuals. After identifying these parameters, the group needed to develop a method for optimizing these data. The group began by using the SOR algorithm to solve a large, sparse system of equations and were able to successfully identify an optimal parameter for the relaxation parameter,\n", " $\\omega$, mathematically and then proceeded to find the optimal value by fixing a particular tolerance level and varying the relaxation parameter in order to see which $\\omega$ value minimized the number of iterations. Graphs of the residuals were generated and the group was able to visualize how varying the relaxation parameter changes the associated residual. The group employed a similar method to find the optimal strength of connection parameter (which cannot be easily found mathematically) and proceeded to investigate how other parameters influence the strength of connection parameter as well as the residual and computation times. Specifically, the group investigated how the strength of connection parameter influenced the post and pre-iterations. The group's approach yielded unclear results regarding the optimal settings for the GMRES solver and different methods should be investigated in the future. The group also computed the optimal parameters by minimizing a function which takes in a desired parameter value (either strength of connection or relaxation) and returns the number of iterations needed to reach the desired tolerance. This was done in an attempt to accelerate and increase the accuracy of the group's optimization and yielded somewhat unclear results (with the $\\theta$ parameter specifically).\n", "\n", "\n", " There are a few main conclusions to draw from the results of the experiments above. First, it is best to note again that all of these experiments were done using sparse, diagonally dominant matrices, these results might vary if given matrices with different properties.\n", "\n", " The first conclusion deals with the pre- and post-smoother iterations. The group initially was unsure of how the pre- and post-smoother iterations would impact the total number of iterations. According to the group's testing, increasing the pre- and post-smoothers to their maximum values appears to decrease the amount of overall iterations needed to reach a solution; increasing the pre- and post-smoothers to their maximum values (or close to their maximum values) resulted in the fewest iterations. This result was seemingly independent of the other settings and were consistent between all experiments. In general, the group's findings suggest that when the pre- and post- iterations increase, the solver iterations decrease. Unfortunately, an error was identified which could not be corrected in the alloted time. The answer for the optimal values differed if the order of the elements in iteration_experiments changed. This insinuates that the group's method is not locating the correct residual and must be addressed. \n", "\n", " The second conclusion deals with the SOR relaxation parameter ($\\omega$). Specifically, the difference between the optimal $\\omega$ for reducing solver iterations and the theoretical optimal $\\omega$. In all experiments, the tests to find the $\\omega$ that reduces the solver iterations to the fewest possible produced multiple possible $\\omega$ values, the program returns the first one it encounters. This means that, the search for the relaxation parameter was not optimal. Despite this, the unique optimal parameter can be reproduced by decreasing the tolerance level. Generally, the group's approach did not initially follow the standard form of a continuous optimization problem but this was addressed towards the end of the project. Defining an objective function with appropriate bounds (which correspond to the parameter being optimized) and minimizing the function using scipy.optimize.minimize_scalar led to a more rigorous and accurate optimization of the relaxation parameter (SOR), but yielded unclear results when applied to the strength of connection parameter (GMRES). Specifically, the value of the optimal parameter would always tend towards the upper bound of the given range which raises doubts regarding the accuracy of this value.\n", "\n", " In regards to the optimal $\\theta$ value, these data do not seem overly-dependent on the $\\theta$ value. Varying the $\\theta$ parameter did not result in large variations in the number of iterations. This is likely due to the type of matrices the solver is using (sparse, diagonally dominant matrices). The data generated by the group suggests that almost all the values of $\\theta$ seem to result in a solution being found in the same number of iterations. It is important to note that varying the post and pre-iteration settings does lead to variations in the optimal strength of connection parameter.\n", "\n", " It is critical to address the discrepancy between the optimal theta calculated using the \"brute-force\" approach (looping over various post/pre-iteration values and $\\theta$) and the optimal $\\theta$ value generated by minimizing an objective function which takes in a value for $\\theta$ and returns the iterations necessary to reach the desired residual. While the latter approach is much quicker, the tendency for the optimal $\\theta$ to consistently be the upper bound is concerning and may have occured for several reasons including: a limitation in the PyAmg source code when building the solver, a quirk in the minimize_scalar() which leads to problems when computing the optimal $\\theta$. Perhaps the optimal parameter actually tends to be the upper bound, however the group couldn't identify a reason for this to be the case. The fact that this approach worked seamlessly for the relaxation parameter (in the SOR solver) suggests that this problem is specific to $\\theta$ which may insinuate an issue in the PyAmg source code and needs to be investigated in the future.\n", "\n", " Towards the very end of the semester, it was discovered that the group's program does not locate the correct residual and iteration count (the ones associated with the optimal parameter) effectively. Essentially, the students' method for locating the aforementioned residual would return the first residual in a set of equal-length residuals. This first residual is generally not the one associated with the optimal parameter. Due to time constraints, the group was unable to rework the necessary code. Unfortunately, theta cannot readily be computed mathematically (as is the case for $\\omega$) so the group is uncertain of the accuracy of their calculated optimal $\\theta$ parameter. If future research is conducted, changing the method described about and comparing the value produced with the optimal $\\theta$ parameter calculated using scipy.optimize.minimize_scalar() may prove to be a useful starting point.\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "9hqHBivLOXgD" }, "source": [ "##5 Continued Research: Machine Learning\r\n", "\r\n", "The next step in AMG optimization which couldn't be implemented this semester comes from the field of AI and machine learning. Much of our work here is formulaic enough that a machine can be trained to find the optimal parameters for the optimal smoothing methods with enough training data. Using the solutions from either NREL's data or our toy data, an algorithm can close in on the optimal $\\theta$ or $\\omega$, removing the need for the time-consuming brute force or exact method.\r\n", "\r\n", "A method for this [supervised learning](https://en.wikipedia.org/wiki/Supervised_learning) would resemble a standard linear regression model:\r\n", "1. Define the input matrix the model will iterate on, the pre-calculated output matrix it will train towards, and whatever parameters the model will vary to get from the input to the output.\r\n", "2. Have the model run our functions while varying the parameters, narrowing the error between the input and output data with each iteration\r\n", "3. When the model has passed a certain error threshold, output the optimal parameters it found.\r\n", "4. Once the model is trained, it can do this without needing to be given the output matrix beforehand.\r\n", "\r\n", "[Keras](https://keras.io/about/) for TensorFlow and [scikit-learn](https://scikit-learn.org/stable/supervised_learning.html#supervised-learning) are the two most accessible frameworks for implementing machine learning in Colab/Jupyter. A problem we ran into was getting a TensorFlow instance to accept both the data matrix and the solution matrix to train against. The former is a Scipy matrix and the latter is a Numpy array, and they will both need to be converted into [tensors](https://colab.research.google.com/github/tensorflow/docs/blob/master/site/en/guide/tensor.ipynb#scrollTo=fDmFtFM7k0R2) (or the Scikit equivalent) for the machine learning instance to operate on them.\r\n", "\r\n", "A subset of machine learning is *deep learning*, which uses multiple inference layers to reach its desired outcome, mimicking a human brain. Each additional layer requires lots of computing power, so while machine learning has been around for decades, deep learning is relatively new. This further accelerates the task by using a [neural network](https://towardsdatascience.com/understanding-neural-networks-from-neuron-to-rnn-cnn-and-deep-learning-cd88e90e0a90) to have the model learn either supervised or unsupervised. Several kinds of neural networks have come about, including artificial ([ANN](https://en.wikipedia.org/wiki/Artificial_neural_network)), convolutional ([CNN](https://en.wikipedia.org/wiki/Convolutional_neural_network)), recurrent ([RNN](https://en.wikipedia.org/wiki/Recurrent_neural_network)), recursive ([RvNN](https://en.wikipedia.org/wiki/Recursive_neural_network)), and the more novel graph ([GNN](https://towardsdatascience.com/a-gentle-introduction-to-graph-neural-network-basics-deepwalk-and-graphsage-db5d540d50b3)). A team at the Weizmann Institute in Israel has [recently researched](https://arxiv.org/pdf/2003.05744.pdf) using a single GNN to find the optimal prolongation operator for any AMG. Their approach differs from ours by using [unsupervised learning](https://en.wikipedia.org/wiki/Unsupervised_learning) to accelerate the prolongation/interpolation step, while ours focuses accelerating the smoothing and iterative method, and would have involved supervised learning for those parameters. Although these approaches are very different, they both accomplish the same goal of solving large lineary systems quickly. Future research may involve replicating or improving this approach, or using it in conjunction with our approach." ] } ] }