optimization_svm.ipynb

{
  "cells": [
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Problema de otimização: SVM"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Implementação baseada no exercício da disciplina de redes neurais. \n",
        "Grupo:\n",
        "- Jorge \n",
        "- Livia\n",
        "- Lucas Cilento\n",
        "- Paulo Viana\n",
        "- Thiago Melo"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Implementação baseada na questão 21.4 do livro \"Neural Networks and Statistical Learning\" dos autores Ke-Lin Du , M. N. S. Swamy. \n",
        "Referências extras:\n",
        "- Pattern Recognition and Machine Learning, Christopher M. Bishop\n",
        "- Method of Lagrange Multipliers: The Theory Behind Support Vector Machines (Tutorial em 3 partes): https://machinelearningmastery.com/method-of-lagrange-multipliers-the-theory-behind-support-vector-machines-part-1-the-separable-case\n",
        "\n",
        "\n",
        "\n",
        "O enunciado do livro para a questão é:"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "ewBiBjs4kkyu"
      },
      "source": [
        "# Livro Questão 21.4\n",
        "For Example (20.1), establish the QP problem, given by (21.7)-(21.9), for SVM learning.\n",
        "Solve $\\alpha_p$, for $p = 1, ..., 4.$"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "O Exemplo 20.1, junto com a função Kernel a ser utilizada:"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "XRRF9s57k_gT"
      },
      "source": [
        "# Livro:  Example (20.1)\n",
        "Define the kernel $k(\\mathbf{x},\\mathbf{x_i}) = (1 + \\mathbf{x^T}\\mathbf{x_i})^2,$ where $\\mathbf{x} = (x_1,x_2)^T$ and $\\mathbf{x_i} = (x_{i1}, x_{i2})^T.$\n",
        "The training samples $\\mathbf{x_1} = (-1,-1)$ and $\\mathbf{x_4} = (+1,+1)$ belong to class 0, and $\\mathbf{x_2} = (-1,+1), \\mathbf{x_3} = (+1,-1)$ to class 1. \n",
        "\n",
        "Expanding the kernel function, we have \n",
        "\n",
        "$k(\\mathbf{x}, \\mathbf{x_i}) = 1 + x_1^2x_{i1}^2 + 2x_1x_2x_{i1}x_{i2} + x_2^2x_{i2}^2 +2x_1x_{i1} +2x_2x_{i2} = \\phi(\\mathbf{x}) \\cdot \\phi(\\mathbf{x_i}),$ \n",
        "\n",
        "where \n",
        "\n",
        "$\\phi(\\mathbf{x}) = (1,x_1^2, \\sqrt{2}x_1x_2, x_2^2, \\sqrt{2}x_1, \\sqrt{2}x_2)^T, $ \n",
        "\n",
        "$\\phi(\\mathbf{x_i}) =  (1, x_{i1}^2, \\sqrt{2}x_{i1}x_{i2}, x_{i2}^2, \\sqrt{2}x_{i1}, \\sqrt{2}x_{i2})^T.$\n",
        "\n",
        "The feature space defined by $\\phi(\\mathbf{x})$ is six-dimensional. To discriminate the four examples in the feature space, we define the decision boundary in $x_1x_2 = 0$. When $x_1x_2 \\ge 0$, an example is categorized into class 0, and otherwise into class 1.\n",
        "\n"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "O Suport Vector Machine Vector também pode ser chamado de maximum margin classifier, por se tratar de um problema de maximização das margens entre a linha que separa linearmente as classes e os pontos que representam os dados. A formulação do problema nos dá um problema de Programação Quadrática, em que temos que minimizar uma função quadrática $f(\\mathbf{w}, \\theta)$, em que $\\mathbf{w}$ são os pesos e $\\theta$ o intercept, obedecendo um conjunto de restrições dadas por:\n",
        "\n",
        "$g_p = y_p(\\mathbf{w}^T \\phi(x_p) + \\theta) \\geq 1$, $p \\in 1,...,N$\n",
        "\n",
        "Em que temos uma restrição desse tipo para cada um dos $N$ pontos.\n",
        "\n",
        "Para resolver o problema, pode-se utilizar o método de multiplicadores de Lagrande, em que definimos uma Lagrangiana:\n",
        "\n",
        "$L (\\mathbf{w}, b, \\mathbf{\\alpha}) = f(\\mathbf{w}, \\theta) - \\sum_p \\alpha_p g_p $\n",
        "\n",
        "Os $\\alpha_p$ são os multiplicadores de Lagrande. Para encontrar a solução ótima, temos que ter as derivadas em relação a $\\mathbf{w}, \\theta$ e $\\alpha$ da Lagrangiana igual a 0. Após aplicar o método e fazendo uma substituição pela função de Kernel $k(\\mathbf{x}, \\mathbf{x_i}) = \\phi(\\mathbf{x}) \\cdot \\phi(\\mathbf{x_i})$, nós obtemos as equações:\n",
        " "
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "c6ZXUsDRrJDY"
      },
      "source": [
        "# Livro: Equações (21.7)-(21.9)\n",
        "\n",
        "- (21.7) $E_{SVM} = \\frac{1}{2} \\sum_{p=1}^N\\sum_{i=1}^N y_py_ik(\\mathbf{x}_p,\\mathbf{x}_i)\\alpha_p\\alpha_i - \\sum_{p=1}^N \\alpha_p$ \n",
        "\n",
        "Restrições:\n",
        "\n",
        "\n",
        "- (21.8) $\\sum_{p=1}^N y_p\\alpha_p = 0 $, \n",
        "- (21.9) $0 \\le \\alpha_p \\le C$, for $p = 1,...,N$, \n",
        "where $\\alpha_p$ is the weight for the kernel corresponding to the p-th example."
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "StKo8lG62oIp"
      },
      "source": [
        "A equação 21.7 é chamada de representação dual e nosso objetivo é minimiza-la.\n",
        "\n",
        "Para a questão temos que encontrar os valores de $α_p$ que minimizem a função $E_{SVM}$ (Definida em 21.7), atendendo às restrições (21.8 )e (21.9). Tudo isso aplicado ao exemplo 20.1. Notemos que as restrições possuem C como um parâmetro livre. Trata-se de um problema de otimização. Para isso, utilizaremos a biblioteca SciPy. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "l6yKAQfXxBVw"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib.colors as pltcolors\n",
        "from scipy import optimize\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qKkgdbW34fVX"
      },
      "source": [
        "Implementando a função de Kernel e colocando os inputs em um vetor $x$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "id": "n7lpp27GUe1k"
      },
      "outputs": [],
      "source": [
        "def kernel(x1, x2):\n",
        "  return (1+np.matmul(x1.T,x2))**2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "cJ-T5isHe0O9"
      },
      "outputs": [],
      "source": [
        "x = np.array([[-1,-1], [-1,1], [1,-1], [1,1]])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fnfq8U344tfe"
      },
      "source": [
        "Podemos imprimir o valor dos kernels para as instâncias"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Cj5AHNzPtIxq",
        "outputId": "69e69e22-5298-4dac-c95f-4f3e84644a06"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[[9. 1. 1. 1.]\n",
            " [1. 9. 1. 1.]\n",
            " [1. 1. 9. 1.]\n",
            " [1. 1. 1. 9.]]\n"
          ]
        }
      ],
      "source": [
        "KernelMatrix = np.zeros((4,4))\n",
        "for i in range(4):\n",
        "  for j in range(4):\n",
        "    KernelMatrix[i][j]=kernel(x[i], x[j])\n",
        "print(KernelMatrix)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t2S5BXLD46Qe"
      },
      "source": [
        "Notamos que o kernel tem uma certa simetria, ele não depende muito dos valores de i e j, dependendo apenas de se eles são iguais ou não. De forma geral $k(x_i, x_j)$ = 9, se $i=j$ ou 1, se $i\\neq j$. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tdeBInMS6BNH"
      },
      "source": [
        "Para a utilização do algoritmo, devemos trabalhar com as classes -1 ou +1, quero dizer que $y_i\\in \\{-1,+1\\}$ Em seguida, os samples 1 e 4, que estão na classe 0 do exemplo 20.1, são codificados para a classe -1 para o nosso problema, e os samples 2 e 3, que eram da classe 1 no exemplo 20.1, ficam na classe +1."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "id": "Rx2PC2whfBKX"
      },
      "outputs": [],
      "source": [
        "y= np.array([-1,1,1,-1])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZcGypYMF7nc9"
      },
      "source": [
        "Agora, vamos implementar o Problema Quadrático, as restrições (21.8) e (21.9), e implementar o vetor gradiente do problema quadrático, em relação a $\\alpha_p$, pois isso ajuda o solver na otimização. Inicialmente fixamos $C=1$, mas também testamos valores maiores, e o resultado do alpha da o mesmo. Quando informamos as restrições, eu também informo o gradiente das restrições"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "id": "vY3piaQRUr--"
      },
      "outputs": [],
      "source": [
        "n=4\n",
        "# Definindo o problema quadrático\n",
        "def QuadraticProblem( alpha):\n",
        "  Error = 0\n",
        "  for p in range(n):\n",
        "    for i in range(n):\n",
        "      Error+= (0.5)*y[p]*y[i]*kernel(x[i], x[p])*alpha[p]*alpha[i]\n",
        "    Error -= alpha[p]\n",
        "  return Error\n",
        "\n",
        "#Derivada em relação a alpha_p do problema quadrático\n",
        "\n",
        "def QuadraticProblemGrad(alpha):\n",
        "  sum = np.zeros(n)\n",
        "\n",
        "  for p in range(n):\n",
        "    for i in range(n):\n",
        "      sum[p]+= y[p]*y[i]*kernel(x[i], x[p])*alpha[i]\n",
        "    sum[p] -= 1\n",
        "  return sum \n",
        "\n",
        "#Restrição (21.8)\n",
        "def constraint1(alpha):\n",
        "  sum = 0\n",
        "  for i in range(n):\n",
        "    sum+= y[i]*alpha[i]\n",
        "  return sum\n",
        "\n",
        "#Restrição (21.9), os bounds\n",
        "C=1\n",
        "b = (0, C)\n",
        "bounds1 =(b,b,b,b)\n",
        "\n",
        "#Informamos a lista de restrições. Também informamos o gradiente da restrição em relação a $alpha_p$\n",
        "constraints_list = ({'type': 'eq',   'fun': lambda a:  constraint1(a), 'jac': lambda a: y})\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JRrJsvG6-xTh"
      },
      "source": [
        "Realizando a otimização"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "id": "cwlhlSEveTj1"
      },
      "outputs": [],
      "source": [
        "optRes = optimize.minimize(fun=lambda a: QuadraticProblem(a),\n",
        "                              x0=np.ones(n), \n",
        "                              method='SLSQP', \n",
        "                              jac=lambda a: QuadraticProblemGrad(a), \n",
        "                              constraints= constraints_list, bounds = bounds1)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Lk0wF7h_iiDr",
        "outputId": "2394302a-8325-48ba-d211-5995aeab30ac"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "     fun: -0.25\n",
              "     jac: array([ 2.22044605e-16, -4.44089210e-16,  4.44089210e-16, -4.44089210e-16])\n",
              " message: 'Optimization terminated successfully.'\n",
              "    nfev: 3\n",
              "     nit: 3\n",
              "    njev: 3\n",
              "  status: 0\n",
              " success: True\n",
              "       x: array([0.125, 0.125, 0.125, 0.125])"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "optRes"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9Uejdp34_k5N"
      },
      "source": [
        "Imprimindo os $\\alpha_p$ resultante"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "khAGr-RVh-x5",
        "outputId": "13b6ea27-1121-43c0-db91-c760e0a0e2a8"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[0.125 0.125 0.125 0.125]\n"
          ]
        }
      ],
      "source": [
        "alpha_resultante = optRes.x\n",
        "print(alpha_resultante)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "95wGMKT5_r9Y"
      },
      "source": [
        "Temos então que $α_p=0.125$, para todos os $p=1,2,3,4$. Podemos variar o valor de C para ver se encontramos diferenças. Fazendo C=10 e C=100, temos:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wxRLL_fCALrp"
      },
      "source": [
        "Para C=10"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8ag-N4u_lsoJ",
        "outputId": "b9219f64-caa2-4167-b0da-3c516f31c635"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[0.125 0.125 0.125 0.125]\n"
          ]
        }
      ],
      "source": [
        "C=10\n",
        "b = (0, C)\n",
        "bounds1 =(b,b,b,b)\n",
        "optRes = optimize.minimize(fun=lambda a: QuadraticProblem(a),\n",
        "                              x0=np.ones(n), \n",
        "                              method='SLSQP', \n",
        "                              jac=lambda a: QuadraticProblemGrad(a), \n",
        "                              constraints= constraints_list, bounds = bounds1)\n",
        "alpha_resultante = optRes.x\n",
        "print(alpha_resultante)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7-qx0IYwAKc4"
      },
      "source": [
        "$α_p = 0.125$  para todos os p."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WaBJSHOBAYIY"
      },
      "source": [
        "Para C=100"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "uawHGvp9AaXZ",
        "outputId": "e685f74c-1a5d-4352-e253-a6c73231bfee"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[0.125 0.125 0.125 0.125]\n"
          ]
        }
      ],
      "source": [
        "C=100\n",
        "b = (0, C)\n",
        "bounds1 =(b,b,b,b)\n",
        "optRes = optimize.minimize(fun=lambda a: QuadraticProblem(a),\n",
        "                              x0=np.ones(n), \n",
        "                              method='SLSQP', \n",
        "                              jac=lambda a: QuadraticProblemGrad(a), \n",
        "                              constraints= constraints_list, bounds = bounds1)\n",
        "alpha_resultante = optRes.x\n",
        "print(alpha_resultante)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bBExH8duAfTU"
      },
      "source": [
        "$α_p = 0.125$  para todos os p. Vemos então que aumentar o valor de C não altera nossa solução. Podemos tentar diminuir C para 0.1"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "JnzPtyZbAJoL",
        "outputId": "b2969064-dffa-42c9-987d-2e9b98c96455"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[0.1 0.1 0.1 0.1]\n"
          ]
        }
      ],
      "source": [
        "C=0.1\n",
        "b = (0, C)\n",
        "bounds1 =(b,b,b,b)\n",
        "optRes = optimize.minimize(fun=lambda a: QuadraticProblem(a),\n",
        "                              x0=np.ones(n), \n",
        "                              method='SLSQP', \n",
        "                              jac=lambda a: QuadraticProblemGrad(a), \n",
        "                              constraints= constraints_list, bounds = bounds1)\n",
        "alpha_resultante = optRes.x\n",
        "print(alpha_resultante)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ERFHzja4CNwa"
      },
      "source": [
        "Vemos que o $α$ bate no teto do bound, 0.1, porém podemos ver abaixo que o valor da função é -0.24,enquanto que para os outros valores de C o valor é -0.25. Ou seja, abaixar o valor de C resultou em uma menor minimização da função, o que não é o desejado"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8FPr0nu0CFJ_",
        "outputId": "4287d321-08d2-434c-e8ea-f8173ff12cf3"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "     fun: -0.24\n",
              "     jac: array([-0.2, -0.2, -0.2, -0.2])\n",
              " message: 'Optimization terminated successfully.'\n",
              "    nfev: 1\n",
              "     nit: 1\n",
              "    njev: 1\n",
              "  status: 0\n",
              " success: True\n",
              "       x: array([0.1, 0.1, 0.1, 0.1])"
            ]
          },
          "execution_count": 19,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "optRes"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8BkiIEz_DwjY"
      },
      "source": [
        "Logo, nossa resposta é $α_p = 0.125$  para todos os p. "
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [
        "ewBiBjs4kkyu",
        "XRRF9s57k_gT"
      ],
      "name": "Questão_04_REDES_NEURAIS.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}