Quasi_Newton algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "np.set_printoptions(precision=4)\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Quasi_Newton(x_0, H_0, H_update_mechanism):\n",
    "        \n",
    "    #Defining the necessary lists...\n",
    "    x = []\n",
    "    x_diff = []\n",
    "\n",
    "    g = []\n",
    "    g_diff = []\n",
    "\n",
    "    d = []\n",
    "    alpha = []\n",
    "    H = []\n",
    "    B = []\n",
    "\n",
    "    x.append(x_0)\n",
    "    H.append(H_0)\n",
    "    B.append(H_0)\n",
    "    \n",
    "    g.append(grad(x[0]))\n",
    "    print(\"g[{0}]:\\n\".format(0), g[0])\n",
    "    \n",
    "    k = 0\n",
    "    while(True):\n",
    "        print(\"\\n\\nIteration \", k+1)        \n",
    "        \n",
    "        if grad_equals_zero(g[k]):    #If none of the elements are true, or, if all of the elements are zero\n",
    "            return x[k] \n",
    "        else:\n",
    "            d.append(-np.dot(H[k], g[k]))\n",
    "            print(\"d[{0}]:\\n\".format(k), d[k])\n",
    "        \n",
    "        # We use the formula for alpha as given by Conjugate direction method,\n",
    "        # since in the quadratic case the Quasi-Newton algortihm is also a  Conjugate direction algorithm\n",
    "        alpha.append(-np.dot(g[k].T, d[k]).item()/(np.dot(np.dot(d[k].T, Q), d[k]).item()))\n",
    "        print(\"alpha[{0}]: \".format(k), alpha[k])\n",
    "        \n",
    "        x_diff.append(alpha[k]*d[k])\n",
    "        print(\"x_diff[{0}]:\\n\".format(k), x_diff[k])\n",
    "        \n",
    "        x.append(x[k] + x_diff[k])\n",
    "        print(\"x[{0}]:\\n\".format(k+1), x[k+1])\n",
    "        \n",
    "        g.append(grad(x[k+1]))\n",
    "        print(\"g[{0}]:\\n\".format(k+1), g[k+1])\n",
    "        \n",
    "        if grad_equals_zero(g[-1]):\n",
    "            return x[-1]\n",
    "                \n",
    "        g_diff.append(g[k+1] - g[k])\n",
    "        print(\"g_diff[{0}]:\\n\".format(k), g_diff[k])\n",
    "                \n",
    "        if H_update_mechanism is \"BFGS\":\n",
    "            H.append(H_BFGS(x_diff, g_diff, H, k))\n",
    "            print(\"H[{0}]:\\n\".format(k+1), H[k+1])\n",
    "        \n",
    "            B.append(B_update(x_diff, g_diff, B, k))\n",
    "            print(\"B[{0}]:\\n\".format(k+1), B[k+1])\n",
    "        \n",
    "        elif H_update_mechanism is \"DFP\":\n",
    "            H.append(H_DFP(x_diff, g_diff, H, k))\n",
    "            print(\"H[{0}]:\\n\".format(k+1), H[k+1])\n",
    "        \n",
    "        elif H_update_mechanism is \"R1\":\n",
    "            H.append(H_rankone(x_diff, g_diff, H, k))\n",
    "            print(\"H[{0}]:\\n\".format(k+1), H[k+1])\n",
    "        \n",
    "        k = k+1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computing H[k+1] given H[k]\n",
    "\n",
    "def H_BFGS(x_diff, g_diff, H, k):\n",
    "    \n",
    "    t1 = np.dot(x_diff[k], x_diff[k].T)/np.dot(x_diff[k].T, g_diff[k]).item()\n",
    "    \n",
    "    coeff_t1 = 1 + np.dot(np.dot(g_diff[k].T, H[k]), g_diff[k]).item()/np.dot(g_diff[k].T, x_diff[k]).item()\n",
    "    \n",
    "    t1 *= coeff_t1\n",
    "    \n",
    "    temp = np.dot(np.dot(H[k], g_diff[k]), x_diff[k].T)\n",
    "    t2 = (temp + temp.T)/np.dot(g_diff[k].T, x_diff[k]).item()\n",
    "    \n",
    "    return H[k] + t1-t2\n",
    "\n",
    "\n",
    "\n",
    "def H_DFP(x_diff, g_diff, H, k):\n",
    "        \n",
    "    t1 = np.dot(x_diff[k], x_diff[k].T)/np.dot(x_diff[k].T, g_diff[k]).item()\n",
    "    \n",
    "    temp = np.dot(H[k], g_diff[k])\n",
    "    t2 = np.dot(temp,temp.T)/np.dot(g_diff[k].T, temp).item()\n",
    "    \n",
    "    return H[k] + t1-t2\n",
    "\n",
    "\n",
    "\n",
    "def H_rankone(x_diff, g_diff, H, k):\n",
    "    \n",
    "    t = x_diff[k] - np.dot(H[k], g_diff[k])\n",
    "    \n",
    "    Nr = np.dot(t, t.T)\n",
    "    Dr = np.dot(g_diff[k].T, t).item()\n",
    "    \n",
    "    return H[k] + Nr/Dr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def B_update(x_diff, g_diff, B, k):\n",
    "    \n",
    "    t1 = np.dot(g_diff[k], g_diff[k].T)/np.dot(g_diff[k].T, x_diff[k]).item()\n",
    "    \n",
    "    temp = np.dot(B[k], x_diff[k])\n",
    "    t2 = np.dot(temp, temp.T)/np.dot(x_diff[k].T, temp).item()\n",
    "    \n",
    "    return B[k] + t1-t2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad(x):\n",
    "    return (np.dot(Q, x) - b)\n",
    "\n",
    "def grad_equals_zero(g):\n",
    "    return not np.any(np.around(g, decimals=4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g[0]:\n",
      " [[ 0]\n",
      " [-1]]\n",
      "\n",
      "\n",
      "Iteration  1\n",
      "d[0]:\n",
      " [[-0.]\n",
      " [ 1.]]\n",
      "alpha[0]:  0.5\n",
      "x_diff[0]:\n",
      " [[-0. ]\n",
      " [ 0.5]]\n",
      "x[1]:\n",
      " [[0. ]\n",
      " [0.5]]\n",
      "g[1]:\n",
      " [[-1.5]\n",
      " [ 0. ]]\n",
      "g_diff[0]:\n",
      " [[-1.5]\n",
      " [ 1. ]]\n",
      "H[1]:\n",
      " [[1.   1.5 ]\n",
      " [1.5  2.75]]\n",
      "B[1]:\n",
      " [[ 5.5 -3. ]\n",
      " [-3.   2. ]]\n",
      "\n",
      "\n",
      "Iteration  2\n",
      "d[1]:\n",
      " [[1.5 ]\n",
      " [2.25]]\n",
      "alpha[1]:  2.0\n",
      "x_diff[1]:\n",
      " [[3. ]\n",
      " [4.5]]\n",
      "x[2]:\n",
      " [[3.]\n",
      " [5.]]\n",
      "g[2]:\n",
      " [[0.]\n",
      " [0.]]\n",
      "\n",
      "The minimizer: x* =\n",
      " [[3.]\n",
      " [5.]]\n"
     ]
    }
   ],
   "source": [
    "#Implement the algorithm...\n",
    "\n",
    "#Eg. 11.4\n",
    "Q = np.array([[5, -3], [-3, 2]])\n",
    "b = np.array([[0, 1]]).T\n",
    "\n",
    "# initial values...\n",
    "x_0 = np.array([[0, 0]]).T\n",
    "H_0 = np.identity(2)\n",
    "\n",
    "print(\"\\nThe minimizer: x* =\\n\", Quasi_Newton(x_0, H_0, \"BFGS\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g[0]:\n",
      " [[2]\n",
      " [2]]\n",
      "\n",
      "\n",
      "Iteration  1\n",
      "d[0]:\n",
      " [[-2.]\n",
      " [-2.]]\n",
      "alpha[0]:  0.6666666666666666\n",
      "x_diff[0]:\n",
      " [[-1.3333]\n",
      " [-1.3333]]\n",
      "x[1]:\n",
      " [[-0.3333]\n",
      " [ 0.6667]]\n",
      "g[1]:\n",
      " [[-0.6667]\n",
      " [ 0.6667]]\n",
      "g_diff[0]:\n",
      " [[-2.6667]\n",
      " [-1.3333]]\n",
      "H[1]:\n",
      " [[0.5 0. ]\n",
      " [0.  1. ]]\n",
      "\n",
      "\n",
      "Iteration  2\n",
      "d[1]:\n",
      " [[ 0.3333]\n",
      " [-0.6667]]\n",
      "alpha[1]:  1.0\n",
      "x_diff[1]:\n",
      " [[ 0.3333]\n",
      " [-0.6667]]\n",
      "x[2]:\n",
      " [[0.]\n",
      " [0.]]\n",
      "g[2]:\n",
      " [[0.]\n",
      " [0.]]\n",
      "\n",
      "The minimizer: x* =\n",
      " [[0.]\n",
      " [0.]]\n"
     ]
    }
   ],
   "source": [
    "#Implement the algorithm...\n",
    "\n",
    "Q = np.array([[2, 0], [0, 1]])\n",
    "b = np.array([[0, 0]]).T\n",
    "\n",
    "# initial values...\n",
    "x_0 = np.array([[1, 2]]).T\n",
    "H_0 = np.identity(2)\n",
    "\n",
    "print(\"\\nThe minimizer: x* =\\n\", Quasi_Newton(x_0, H_0, \"R1\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g[0]:\n",
      " [[1]\n",
      " [2]]\n",
      "\n",
      "\n",
      "Iteration  1\n",
      "d[0]:\n",
      " [[-1.]\n",
      " [-2.]]\n",
      "alpha[0]:  0.3333333333333333\n",
      "x_diff[0]:\n",
      " [[-0.3333]\n",
      " [-0.6667]]\n",
      "x[1]:\n",
      " [[-0.3333]\n",
      " [ 0.3333]]\n",
      "g[1]:\n",
      " [[-0.6667]\n",
      " [ 0.3333]]\n",
      "g_diff[0]:\n",
      " [[-1.6667]\n",
      " [-1.6667]]\n",
      "H[1]:\n",
      " [[ 0.6222 -0.4222]\n",
      " [-0.4222  0.8222]]\n",
      "B[1]:\n",
      " [[2.4667 1.2667]\n",
      " [1.2667 1.8667]]\n",
      "\n",
      "\n",
      "Iteration  2\n",
      "d[1]:\n",
      " [[ 0.5556]\n",
      " [-0.5556]]\n",
      "alpha[1]:  0.6\n",
      "x_diff[1]:\n",
      " [[ 0.3333]\n",
      " [-0.3333]]\n",
      "x[2]:\n",
      " [[0.]\n",
      " [0.]]\n",
      "g[2]:\n",
      " [[0.]\n",
      " [0.]]\n",
      "\n",
      "The minimizer: x* =\n",
      " [[0.]\n",
      " [0.]]\n"
     ]
    }
   ],
   "source": [
    "#Implement the algorithm...\n",
    "\n",
    "Q = np.array([[3, 1], [1, 2]])\n",
    "b = np.array([[0, 0]]).T\n",
    "\n",
    "# initial values...\n",
    "x_0 = np.array([[0, 1]]).T\n",
    "H_0 = np.identity(2)\n",
    "\n",
    "print(\"\\nThe minimizer: x* =\\n\", Quasi_Newton(x_0, H_0, \"BFGS\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}