5-quasinewton.ipynb

{
 "metadata": {
  "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.9.1-final"
  },
  "orig_nbformat": 2,
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.9.1 64-bit",
   "metadata": {
    "interpreter": {
     "hash": "0d008810b9c8467bcb3ca39aa2180e5b81b3a9acb136aab30d47954377cc5120"
    }
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.core.fromnumeric import argmin\n",
    "import functions as func\n",
    "import plots as mypt\n",
    "import time\n",
    "\n",
    "def is_pos_def(x):\n",
    "    return np.all(np.linalg.eigvals(x) > 0)\n",
    "'''\n",
    "def trust_region_subproblem(gradient, hessian, lambda_0, delta, max_iter=1000):\n",
    "    d = hessian.shape[0]\n",
    "    lambda_j, vector = np.linalg.eig(hessian)\n",
    "    lambda_1 = min(lambda_j)\n",
    "    min_index = argmin(lambda_j)\n",
    "    lam = lambda_0\n",
    "    p = np.linalg.inv(hessian) @ (-gradient) #initial p\n",
    "    count = 0\n",
    "    #print(hessian)\n",
    "    while lam >= -lambda_1:\n",
    "        count += 1\n",
    "        B_ = hessian + lam * np.identity(d)\n",
    "        # try:\n",
    "        if vector[:, min_index].T @ gradient == 0:\n",
    "            p = hard_case(gradient, hessian, delta, lambda_j, vector, min_index)\n",
    "            return p\n",
    "        else:\n",
    "            p = np.linalg.inv(B_) @ (-gradient)\n",
    "        R = np.linalg.cholesky(B_)\n",
    "        q = np.linalg.inv(R.T) @ p\n",
    "\n",
    "        lam_new = lam + (np.linalg.norm(p) / np.linalg.norm(q)) ** 2 * (\n",
    "            (np.linalg.norm(p) - delta) / delta\n",
    "        )\n",
    "        #print(lam)\n",
    "        if lam_new == lam or (count >= max_iter and np.linalg.norm(p) <= delta) or (lam_new < -lambda_1):\n",
    "            # return p\n",
    "            break\n",
    "        else:\n",
    "            lam = lam_new\n",
    "    print(count)\n",
    "    return p\n",
    "'''\n",
    "def sub_solve_2(gx, hx, delta):\n",
    "    atol = 1e-8\n",
    "    gx_norm = np.linalg.norm(gx)\n",
    "    lmds, Q = np.linalg.eig(hx)\n",
    "    minIdx = np.argmin(lmds)\n",
    "    #is_psd  = lmds[minIdx] >= -atol\n",
    "    if Q[:, minIdx].T @ gx != 0:\n",
    "        if lmds[minIdx] > atol:\n",
    "            lmd = 0\n",
    "        else:\n",
    "            lmd = -lmds[minIdx] + atol\n",
    "        while True:\n",
    "            #B_inv = Q @ np.diag(1 / (lmds + lmd)) @ Q.T\n",
    "            B_inv = Q @ np.diag(1 / (lmds + lmd + atol)) @ Q.T\n",
    "            p     = B_inv @ -gx\n",
    "            p_n   = np.linalg.norm(p)\n",
    "            if p_n <= delta + atol:\n",
    "                break\n",
    "            p_n2  = p_n*p_n\n",
    "            q_n2  = p.T @ B_inv @ p\n",
    "            lmd   = lmd + (p_n2/q_n2) * (p_n - delta) / delta\n",
    "    else:\n",
    "        #print(f\"Hard case! (||gx|| = {gx_norm})\")\n",
    "        p_n = delta\n",
    "        if gx_norm == 0:\n",
    "            p = delta*Q[:, minIdx]\n",
    "        else:\n",
    "            p = delta*(-gx/gx_norm)\n",
    "            gBg = gx.T @ hx @ gx\n",
    "            if gBg > 0:\n",
    "                tau = (gx_norm*gx_norm*gx_norm) / (delta * gBg)\n",
    "                if tau < 1:\n",
    "                    p, p_n = p*tau, p_n*tau\n",
    "    return p\n",
    "'''\n",
    "def hard_case(gradient, hessian, delta, lambda_j, vector, min_index):\n",
    "    temp = gradient.T @ hessian @ gradient\n",
    "    tau = min(1, np.linalg.norm(gradient) ** 3 / (delta * temp))\n",
    "    z = vector[:, min_index]  # eigenvector of B corresponding to lamda_1\n",
    "    lamda_1 = lambda_j[min_index]\n",
    "    p = 0\n",
    "    for i in range(np.shape(lambda_j)[0]):\n",
    "        if lamda_1 != lambda_j[i]:\n",
    "            p += -((vector[:, i].T @ gradient) / (lambda_j[i] - lamda_1) * vector[:, i])\n",
    "    p += tau * z\n",
    "    return p\n",
    "\n",
    "def cauchy_point(gradient, hessian, delta):\n",
    "    temp = gradient.T @ hessian @ gradient\n",
    "    if temp <= 0:\n",
    "        tau = 1\n",
    "    else:\n",
    "        tau = min(1, np.linalg.norm(gradient) ** 3 / (delta * temp))\n",
    "    return -tau * delta / np.linalg.norm(gradient) * gradient\n",
    "'''\n",
    "def trust_region(_x, fun, fun_g, fun_h, _delta, delta_max, eta=1e-4, max_iter=1000, tol=1e-6):\n",
    "\txs = [_x] #store x_0\n",
    "\tdeltas = [_delta] #store delta_0\n",
    "\t#Bs = [np.identity(len(_x))] #store B_0\n",
    "\tx = _x #initialize x_k to x_0\n",
    "\tdelta = _delta #initialize delta_k to delta_0\n",
    "\tB = np.identity(len(_x)) #initialize B_k as identity matrix B_0\n",
    "\tcount = 0 #store iterations number\n",
    "\twhile np.linalg.norm(fun_g(x))>tol: #until reach the tolerance\n",
    "\t\tif(count >= max_iter): break #until reach the maximum iterations\n",
    "\t\tcount += 1\n",
    "\t\tg = fun_g(x)\n",
    "\t\t#B = fun_h(x) #this is same as last time in assignment 4\n",
    "\t\ts = sub_solve_2(g, B, delta)\n",
    "\t\t'''\n",
    "\t\tif is_pos_def(B):\n",
    "\t\t\tif np.linalg.norm(np.linalg.inv(B) @ g) <= delta:\n",
    "\t\t\t\ts = -np.linalg.inv(B) @ g\n",
    "\t\t\telse:\n",
    "\t\t\t\ts = trust_region_subproblem(g, B, 0, delta)\n",
    "\t\telse:\n",
    "\t\t\ts = trust_region_subproblem(g, B, -min(np.linalg.eigvals(B)+1e-8), delta)  #set the first lambda to almost -lambda_1\n",
    "\t\t'''\n",
    "\t\ty = fun_g(x+s) - g\n",
    "\t\tared = fun(x) - fun(x+s) #actual reduction\n",
    "\t\tpred = -(g.T @ s + (s.T @ B @ s)/2) #predicted reduction\n",
    "\t\trho = ared / pred\n",
    "\t\t#print(\"ared:\",ared,\"\\tpred:\",pred,\"\\tdelta:\",delta,\"\\t||g||\",np.linalg.norm(fun_g(x)))\n",
    "\t\tif rho > eta:\n",
    "\t\t\tx = x + s\n",
    "\t\telse:\n",
    "\t\t\tpass #x_k+1=x_k\n",
    "\t\tif rho > 0.75:\n",
    "\t\t\tif np.linalg.norm(s) <= 0.8*delta:\n",
    "\t\t\t\tpass\n",
    "\t\t\telse:\n",
    "\t\t\t\tdelta=2*delta\n",
    "\t\telif 0.1 <= rho <= 0.75:\n",
    "\t\t\tpass\n",
    "\t\telse:\n",
    "\t\t\tdelta=0.5*delta\n",
    "\t\tr=1e-8; tmp=y-B@s\n",
    "\t\tif np.array_equal(y, B @ s):\n",
    "\t\t\tpass\n",
    "\t\telif np.abs(s.T @ tmp) < 1e-8:\n",
    "\t\t\tpass\n",
    "\t\telif np.abs(s.T @ tmp)>=r*np.linalg.norm(s)*np.linalg.norm(tmp): #6.26\n",
    "\t\t\tB=B+np.outer(tmp,tmp)/(tmp.T @ s) #6.24, SR1\n",
    "\t\telse:\n",
    "\t\t\tpass\n",
    "\t\txs.append(x)\n",
    "\t\tdeltas.append(delta)\n",
    "\treturn np.array(xs), fun(x), count, np.array(deltas)\n",
    "\n",
    "def test(fun, fun_g, fun_h, dim=2, name=\"\",tol=1e-6):\n",
    "    time_s,time_f = 0,0\n",
    "    iter_s,iter_f = 0,0\n",
    "    succ_s = 0\n",
    "    val_s,val_f = 0,0\n",
    "    x_s,x_f = 0,0\n",
    "    points = 100\n",
    "    for i in range(points):\n",
    "        #x = np.random.rand(length) * 5\n",
    "        x = np.zeros(dim)\n",
    "        for j in range(dim):\n",
    "            x[j] = np.random.uniform(-5, 5)\n",
    "        tmp1 = time.process_time()\n",
    "        xs, val, c, dels = trust_region(x, fun, fun_g, fun_h,0.1,10,tol=tol)\n",
    "        tmp2 = time.process_time()\n",
    "        x_ = xs[-1]\n",
    "        if np.linalg.norm(fun_g(x_)) <= tol:\n",
    "            succ_s += 1\n",
    "            time_s += tmp2-tmp1\n",
    "            iter_s += c\n",
    "            #x_s += x_\n",
    "            val_s += val\n",
    "        else:\n",
    "            time_f += tmp2-tmp1\n",
    "            iter_f += c\n",
    "            #x_f += x_\n",
    "            val_f += val\n",
    "    succ_f=points-succ_s\n",
    "    '''\n",
    "    if(succ_s>0 and succ_f>0):\n",
    "        print(name+\" succ\", \"\\ttime:\", round(time_s/succ_s,6), \"\\titer:\", round(iter_s/succ_s,4), \"\\ty:\", val_s/succ_s, \"\\tx:\", x_s/succ_s, \"\\tsucc rate:\", succ_s/points)\n",
    "        print(name+\" fail\", \"\\ttime:\", round(time_f/succ_f,6), \"\\titer:\", round(iter_f/succ_f,4), \"\\ty:\", val_f/succ_f, \"\\tx:\", x_f/succ_f, \"\\tfail rate:\", succ_f/points)\n",
    "    elif(succ_s==0):\n",
    "        #print(name+\" without successful convergence.\")\n",
    "        print(name+\" fail\", \"\\ttime:\", round(time_f/succ_f,6), \"\\titer:\", round(iter_f/succ_f,4), \"\\ty:\", val_f/succ_f, \"\\tx:\", x_f/succ_f, \"\\tfail rate:\", succ_f/points)\n",
    "    elif(succ_f==0):\n",
    "        print(name+\" succ\", \"\\ttime:\", round(time_s/succ_s,6), \"\\titer:\", round(iter_s/succ_s,4), \"\\ty:\", val_s/succ_s, \"\\tx:\", x_s/succ_s, \"\\tsucc rate:\", succ_s/points)\n",
    "        #print(name+\" without failed convergence.\")\n",
    "    '''\n",
    "    if(succ_s>0 and succ_f>0):\n",
    "        print(\"& succ\", \"&\", round(time_s/succ_s,6), \"&\", round(iter_s/succ_s,4), \"&\", succ_s/points)\n",
    "        print(\"& fail\", \"&\", round(time_f/succ_f,6), \"&\", round(iter_f/succ_f,4), \"&\", succ_f/points)\n",
    "    elif(succ_s==0):\n",
    "        #print(name+\" without successful convergence.\")\n",
    "        print(\"& fail\", \"&\", round(time_f/succ_f,6), \"&\", round(iter_f/succ_f,4), \"&\", succ_f/points)\n",
    "    elif(succ_f==0):\n",
    "        print(\"& succ\", \"&\", round(time_s/succ_s,6), \"&\", round(iter_s/succ_s,4), \"&\", succ_s/points)\n",
    "        #print(name+\" without failed convergence.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "'\\nx0 = np.array([-1.8, 2.1])\\n\\nX1,X2,Y=mypt.gen2d(0.01,-3.5,3.5,2) #generate points set for the contour plot\\n\\nf1_xs, f1_y, f1_count, f1_dels = trust_region(np.array(x0), func.f1, func.g1, func.h1, 0.1, 10)\\nf2_xs, f2_y, f2_count, f2_dels = trust_region(np.array(x0), func.f2, func.g2, func.h2, 0.1, 10)\\nf3_xs, f3_y, f3_count, f3_dels = trust_region(np.array(x0), func.f3, func.g3, func.h3, 0.1, 10)\\nf4_xs, f4_y, f4_count, f4_dels = trust_region(np.array(x0), func.f4, func.g4, func.h4, 0.1, 10)\\nf5_xs, f5_y, f5_count, f5_dels = trust_region(np.array(x0), func.f4, func.g4, func.h4, 0.1, 10)\\n\\nmypt.plot_deltas(f1_dels, name=\"f1\")\\nmypt.plot_deltas(f2_dels, name=\"f2\")\\nmypt.plot_deltas(f3_dels, name=\"f3\")\\nmypt.plot_deltas(f4_dels, name=\"f4\")\\nmypt.plot_deltas(f5_dels, name=\"f5\")\\n'"
      ]
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "source": [
    "'''\n",
    "x0 = np.array([-1.8, 2.1])\n",
    "\n",
    "X1,X2,Y=mypt.gen2d(0.01,-3.5,3.5,2) #generate points set for the contour plot\n",
    "\n",
    "f1_xs, f1_y, f1_count, f1_dels = trust_region(np.array(x0), func.f1, func.g1, func.h1, 0.1, 10)\n",
    "f2_xs, f2_y, f2_count, f2_dels = trust_region(np.array(x0), func.f2, func.g2, func.h2, 0.1, 10)\n",
    "f3_xs, f3_y, f3_count, f3_dels = trust_region(np.array(x0), func.f3, func.g3, func.h3, 0.1, 10)\n",
    "f4_xs, f4_y, f4_count, f4_dels = trust_region(np.array(x0), func.f4, func.g4, func.h4, 0.1, 10)\n",
    "f5_xs, f5_y, f5_count, f5_dels = trust_region(np.array(x0), func.f4, func.g4, func.h4, 0.1, 10)\n",
    "\n",
    "mypt.plot_deltas(f1_dels, name=\"f1\")\n",
    "mypt.plot_deltas(f2_dels, name=\"f2\")\n",
    "mypt.plot_deltas(f3_dels, name=\"f3\")\n",
    "mypt.plot_deltas(f4_dels, name=\"f4\")\n",
    "mypt.plot_deltas(f5_dels, name=\"f5\")\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "'\\nmypt.plot_convergence(f1_xs,func.f1,name=\"f1\")\\nmypt.plot_convergence(f2_xs,func.f2,name=\"f2\")\\nmypt.plot_convergence(f3_xs,func.f3,name=\"f3\")\\nmypt.plot_convergence(f4_xs,func.f4,name=\"f4\")\\nmypt.plot_convergence(f5_xs,func.f5,name=\"f5\")\\n'"
      ]
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "source": [
    "'''\n",
    "mypt.plot_convergence(f1_xs,func.f1,name=\"f1\")\n",
    "mypt.plot_convergence(f2_xs,func.f2,name=\"f2\")\n",
    "mypt.plot_convergence(f3_xs,func.f3,name=\"f3\")\n",
    "mypt.plot_convergence(f4_xs,func.f4,name=\"f4\")\n",
    "mypt.plot_convergence(f5_xs,func.f5,name=\"f5\")\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "'\\nmypt.plot2d(X1,X2,Y,x0,func.f1,f1_xs,f1_count,f1_y,name=\"f1\")\\nmypt.plot2d(X1,X2,Y,x0,func.f2,f2_xs,f2_count,f2_y,name=\"f2\")\\nmypt.plot2d(X1,X2,Y,x0,func.f3,f3_xs,f3_count,f3_y,name=\"f3\")\\nmypt.plot2d(X1,X2,Y,x0,func.f4,f4_xs,f4_count,f4_y,name=\"f4\")\\nmypt.plot2d(X1,X2,Y,x0,func.f5,f5_xs,f5_count,f5_y,name=\"f5\")\\n'"
      ]
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "source": [
    "'''\n",
    "mypt.plot2d(X1,X2,Y,x0,func.f1,f1_xs,f1_count,f1_y,name=\"f1\")\n",
    "mypt.plot2d(X1,X2,Y,x0,func.f2,f2_xs,f2_count,f2_y,name=\"f2\")\n",
    "mypt.plot2d(X1,X2,Y,x0,func.f3,f3_xs,f3_count,f3_y,name=\"f3\")\n",
    "mypt.plot2d(X1,X2,Y,x0,func.f4,f4_xs,f4_count,f4_y,name=\"f4\")\n",
    "mypt.plot2d(X1,X2,Y,x0,func.f5,f5_xs,f5_count,f5_y,name=\"f5\")\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "'\\ntest(func.f1, func.g1, func.h1, name=\"f1\")\\ntest(func.f2, func.g2, func.h2, name=\"f2\")\\ntest(func.f3, func.g3, func.h3, name=\"f3\")\\ntest(func.f4, func.g4, func.h4, name=\"f4\")\\ntest(func.f5, func.g5, func.h5, name=\"f5\")\\n'"
      ]
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "source": [
    "'''\n",
    "test(func.f1, func.g1, func.h1, name=\"f1\")\n",
    "test(func.f2, func.g2, func.h2, name=\"f2\")\n",
    "test(func.f3, func.g3, func.h3, name=\"f3\")\n",
    "test(func.f4, func.g4, func.h4, name=\"f4\")\n",
    "test(func.f5, func.g5, func.h5, name=\"f5\")\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "( 1.0 1.0 1.0 1.0 ) & succ & 0.001563 & 6.35 & 1.0\n",
      "( 1.0 2.0 65.0 730.0 ) & succ & 0.002656 & 6.99 & 1.0\n",
      "( 1.0 4.0 8257.0 1063612.0 ) & succ & 0.002344 & 7.98 & 1.0\n",
      "( 1.0 7.0 794689.0 1163325079.0 ) & succ & 0.002656 & 8.24 & 1.0\n",
      "( 1.0 11.0 67903553.0 1130881471003.0 ) & succ & 0.002656 & 9.21 & 1.0\n",
      "( 1.0 16.0 5436612673.0 1030586541944248.0 ) & succ & 0.002188 & 9.78 & 1.0\n",
      "( 1.0 22.0 417753473089.0 9.01598398323939e+17 ) & succ & 0.004577 & 25.8485 & 0.99\n",
      "& fail & 0.125 & 1000.0 & 0.01\n"
     ]
    }
   ],
   "source": [
    "x_test = np.array([1,-1.6,0.9,1.2]) #4d starting point\n",
    "dim = x_test.shape[0]\n",
    "x_min = np.zeros(dim) #4d minimum point, with result is 0\n",
    "#print(x_test,x_min)\n",
    "\n",
    "Shape = np.identity(dim)\n",
    "\n",
    "for i in range(7):\n",
    "    \n",
    "    for j in range(dim):\n",
    "        for k in range(dim):\n",
    "            if j==k:\n",
    "                Shape[j,k] = Shape[j,k]+i*(j*k)**(3*i) #Shape transforms by changing eigenvalues\n",
    "    \n",
    "    #dim=i\n",
    "    #Shape = np.identity(dim)\n",
    "    #print(Shape)\n",
    "    #Shape[dim-1,dim-1] = 10\n",
    "    #print(Shape)\n",
    "    print('(',np.linalg.eig(Shape)[0][0],np.linalg.eig(Shape)[0][1],np.linalg.eig(Shape)[0][2],np.linalg.eig(Shape)[0][3],') ',end='')\n",
    "    f6 = lambda x: x.T @ Shape @ x #ellipsoid function, min when all x=0\n",
    "    g6 = lambda x: Shape @ x + Shape.T @ x\n",
    "    h6 = 0 #don't need\n",
    "    #print('minimum:',f6(x_min))\n",
    "    #print(dim,'',end='')\n",
    "    test(f6, g6, h6, dim, name=\"f6\") #test on different dimensions\n",
    "    '''\n",
    "    for j in range(100):\n",
    "        tmp1 = time.process_time()\n",
    "        xs, val, c, dels = trust_region(x_test, f6, g6, h6,0.1,10,tol=1e-6)\n",
    "        tmp2 = time.process_time()\n",
    "        x_ = xs[-1]\n",
    "    #print(i+1,end='')\n",
    "    if np.linalg.norm(g6(x_)) <= 1e-6:\n",
    "        print(\"\\tsucc.\\t\",end='')\n",
    "    print('minimize:',val,'\\ttime:',tmp2-tmp1,'\\titerations:',c)\n",
    "    '''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ]
}