Skip to content

Instantly share code, notes, and snippets.

@tuf22191

tuf22191/Conjugate Gradient Method 6 20 2018

Created June 20, 2018 09:52
Show Gist options
  • Save tuf22191/5c48279b364a9ed03a28a0566222eb93 to your computer and use it in GitHub Desktop.
Save tuf22191/5c48279b364a9ed03a28a0566222eb93 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
Conjugate Gradient Method 6 20 2018
{
"cells": [
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Today we will be talking about the conjugate gradient method. This method is a great algorithm/way to find the optimal solutions to quadratic programming problems. The standard form of the quadratic programming problem is this (as seen https://en.wikipedia.org/wiki/Quadratic_programming): \n",
"\n",
" \n",
"\\begin{array}{ll}\n",
"\\text{minimize} & \\frac{1}{2}x^T Ax + c^Tx \\\\\n",
"\\text{subject to}& A x \\le b \\\\\n",
"\\end{array}\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As the article discusses, there are many ways to solve this problem, but we will be using the conjugate gradient method. A good explanation of the algorithm here: \n",
"https://www.youtube.com/watch?v=lI4pm6VGA-Y\n",
"\n",
"Anyways, without getting bogged down into the details of why it works. We will for the time being only care about the \"how\" part - as in, how do we program this algorithm? In future articles we might add on and talk about why the algorithm works. \n",
"\n",
"The below algorithms is taken directly from Wikipedia: https://en.wikipedia.org/wiki/Conjugate_gradient_method\n",
"\n",
"\\begin{align}\n",
"& \\mathbf{r}_0 := \\mathbf{b} - \\mathbf{A x}_0 \\\\\n",
"& \\mathbf{p}_0 := \\mathbf{r}_0 \\\\\n",
"& k := 0 \\\\\n",
"& \\text{repeat} \\\\\n",
"& \\qquad \\alpha_k := \\frac{\\mathbf{r}_k^\\mathsf{T} \\mathbf{r}_k}{\\mathbf{p}_k^\\mathsf{T} \\mathbf{A p}_k} \\\\\n",
"& \\qquad \\mathbf{x}_{k+1} := \\mathbf{x}_k + \\alpha_k \\mathbf{p}_k \\\\\n",
"& \\qquad \\mathbf{r}_{k+1} := \\mathbf{r}_k - \\alpha_k \\mathbf{A p}_k \\\\\n",
"& \\qquad \\hbox{if } r_{k+1} \\text{ is sufficiently small, then exit loop} \\\\\n",
"& \\qquad \\beta_k := \\frac{\\mathbf{r}_{k+1}^\\mathsf{T} \\mathbf{r}_{k+1}}{\\mathbf{r}_k^\\mathsf{T} \\mathbf{r}_k} \\\\\n",
"& \\qquad \\mathbf{p}_{k+1} := \\mathbf{r}_{k+1} + \\beta_k \\mathbf{p}_k \\\\\n",
"& \\qquad k := k + 1 \\\\\n",
"& \\text{end repeat} \\\\\n",
"& \\text{The result is } \\mathbf{x}_{k+1}\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#We will use the iterative approach of the algorithm, since it scales well \n",
"#however this method is indeed APPROXIMATE\n",
"\n",
"# A_mat = np.array([[2 , -1, 0],[-1, 2, -1],[0,-1,2]]) # a positive semidefinate array\n",
"# dimA = len(A_mat)\n",
"# b_vect = np.array([1,2,3])\n",
"# b_vect = b_vect.reshape((dimA,1))\n",
"\n",
"\n",
"def congugateGradient_Iterative_AND_Approximate(A_mat, b_vect):\n",
"#key to remember is that x_n is a vector that is Nx1\n",
"#remember that the only requirement is that A be a positive definate matrix \n",
" dimA = len(A_mat)\n",
" x0 = np.array([np.random.sample()*100 for x in range(0, len(A_mat))])\n",
" #x0 = np.array([50.,20.,100.])\n",
" # x0 = np.array([2.,1.])\n",
"\n",
" x0 = x0.reshape((dimA,1))\n",
" #print(str(x0) + \" is x0\")\n",
" #r0, k0, p0\n",
" r_k = b_vect - np.dot(A_mat,x0).reshape((dimA, 1))\n",
" #print(str(r_k) + ' is r0')\n",
" k=0\n",
" p_k = np.array(r_k)\n",
" #print(str(p_k) + \" is p0\")\n",
" x_k = np.array(x0)\n",
"\n",
" rvalues = [r_k]\n",
" pvalues = [p_k]\n",
" xvalues = [x_k]\n",
" alpha_values = []\n",
" bvalues = []\n",
" gradient_steps = [] # where we store the intermediate steps to the optimum value\n",
" tots_max = dimA+1 # incase our program misbehaves (there is a bug), we want to end the while loop\n",
" count = 0\n",
"\n",
" while(True): #while it is not sufficiently small\n",
" r_k = rvalues[k]\n",
" p_k = pvalues[k]\n",
" x_k = xvalues[k]\n",
"\n",
" alpha_k = np.dot(r_k.T, r_k)/(np.dot(p_k.T, np.dot(A_mat, p_k)))\n",
" alpha_k = alpha_k[0][0]\n",
" x_kplus1 = x_k+alpha_k*p_k\n",
" r_kplus1 = r_k-alpha_k*np.dot(A_mat, p_k)\n",
"\n",
"\n",
" if(np.sqrt(np.dot(r_kplus1.T, r_kplus1)[0][0])<1e-30):\n",
" print(str(k+1) + ' steps to break')\n",
" alpha_values.append(alpha_k)\n",
" break\n",
"\n",
" beta_k = np.dot(r_kplus1.T, r_kplus1)/np.dot(r_k.T, r_k)\n",
" beta_k =beta_k[0][0]\n",
" p_kplus1 = r_kplus1 + beta_k*p_k\n",
"\n",
" #update the iterations\n",
" k = k+1\n",
" rvalues.append(r_kplus1)\n",
" xvalues.append(x_kplus1)\n",
" pvalues.append(p_kplus1)\n",
" alpha_values.append(alpha_k)\n",
" bvalues.append(beta_k)\n",
" count+=1\n",
" if(count >=tots_max):\n",
" break\n",
" \n",
" return xvalues, pvalues, alpha_values, rvalues, bvalues "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we are going to use the exact example from the Wikipedia article! Go here: https://en.wikipedia.org/wiki/Conjugate_gradient_method and scroll to the \"numerical example\". Let's see if we can get the solution of [0.09090909, .63636364] transposed for the answer. "
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The approximate solution to the problem is: \n",
"[[ 0.09090909]\n",
" [ 0.63636364]]\n"
]
}
],
"source": [
"# print(rvalues)\n",
"# print(pvalues)\n",
"# print(xvalues)\n",
"A_mat = np.array([[4. , 1.],[1., 3.]]) # a positive semidefinate array\n",
"b_vect = np.array([1.,2.])\n",
"b_vect = b_vect.reshape((len(A_mat),1))\n",
"xvalues, pvalues, alpha_values, rvalues, bvalues = congugateGradient_Iterative_AND_Approximate(A_mat, b_vect)\n",
"\n",
"debug = False \n",
"if(debug == True):\n",
" #for troubleshooting the values \n",
" for index, x_ in enumerate(xvalues):\n",
" print(str(x_) + \" - \" + 'x'+str(index))\n",
" for index, x_ in enumerate(pvalues):\n",
" print(str(x_) + \" - \" + 'p'+str(index))\n",
" for index, x_ in enumerate(alpha_values):\n",
" print(str(x_) + \" - \" + 'a'+str(index))\n",
" for index, x_ in enumerate(rvalues):\n",
" print(str(x_) + \" - \" + 'r'+str(index)) \n",
" for index, x_ in enumerate(bvalues):\n",
" print(str(x_) + \" - \" + 'b'+str(index)) \n",
"\n",
"#the answer\n",
"print('The approximate solution to the problem is: ')\n",
"print(xvalues[-1]+alpha_values[-1]*pvalues[-1])\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course this is only an approximate solution! Again although we got the exact numbers the Wikipedia article got, we can visit WolframAlpha here to see the real answer. \n",
"\n",
"https://www.wolframalpha.com/input/?i=minimize++4x%5E2%2B2xy%2B3y%5E2-x-2y\n",
"\n",
"Looking at the answer, we aren't too close, but we still have an approximation!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can use this solver to solve all sorts of quadratic programming problems!\n",
"Like training support vector machines! Also another example is: \n",
"http://www.uio.no/studier/emner/matnat/ifi/INF-MAT4350/h08/undervisningsmateriale/chap15slides.pdf"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The approximate solution to the problem is: \n",
"[[ 0.66666667]\n",
" [ 0.33333333]]\n"
]
}
],
"source": [
"A_mat = np.array([[2. , -1.],[-1., 2.]]) # a positive semidefinate array\n",
"b_vect = np.array([1.,0.])\n",
"b_vect = b_vect.reshape((len(A_mat),1))\n",
"xvalues, pvalues, alpha_values, rvalues, bvalues = congugateGradient_Iterative_AND_Approximate(A_mat, b_vect)\n",
"print('The approximate solution to the problem is: ')\n",
"print(xvalues[-1]+alpha_values[-1]*pvalues[-1])\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.14"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment