qr.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Decomposição QR e Algoritmo QR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O objetivo deste notebook é ilustrar e explicar brevemente a minha implementação do método de decomposição QR,\n",
    "usando o método $\\textit{Gram-Schmidt Orthogonalization}$, e a posterior construção do algoritmo QR, para\n",
    "encontrar os $\\lambda' s$, que são os autovalores da matriz cuja decomposição QR foi encontrada."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos considerar o caso de uma matriz $A$ simétrica.\n",
    "\n",
    "Vamos gerar primeiro uma matriz $M$ cujas entradas são geradas de uma distribuição Gaussiana padronizada,\n",
    "e somaremos $M + M^{T}$, para obtermos a matriz simétrica $A$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "#para fins de reprodutibilidade\n",
    "np.random.seed(98489)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.60276449  0.09178402  0.44795222]\n",
      " [-1.0051296   0.2770554   0.87940038]\n",
      " [-0.86943815 -0.68699472  0.03559767]]\n"
     ]
    }
   ],
   "source": [
    "M = np.random.normal(size=(3,3))\n",
    "print (M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.20552897 -0.91334558 -0.42148593]\n",
      " [-0.91334558  0.5541108   0.19240567]\n",
      " [-0.42148593  0.19240567  0.07119534]]\n"
     ]
    }
   ],
   "source": [
    "A = M + M.T\n",
    "print (A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vou definir agora a minha implementação do método QR, \n",
    "que foi baseada no livro $\\textbf{Numerical Algorithms}$, do $\\textit{Justin Solomon}$ (os comentários foram extraídos do próprio livro:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# numerically unstable\n",
    "def gram_schmidt_qr(A):\n",
    "    m, n = A.shape\n",
    "    Q = np.zeros_like(A)\n",
    "    for i in range(0,n):\n",
    "        # Q\n",
    "        u_i = A[:,i] - ( np.sum( A[:,i].dot(Q[:,:i])*Q[:,:i], axis=1) )\n",
    "        Q[:,i] = u_i/np.linalg.norm(u_i)\n",
    "\n",
    "        # When the columns of A are linearly independent,\n",
    "        # one way to find R is:\n",
    "        # Q.T * A = Q.T * Q * R = R\n",
    "        R = Q.T.dot(A)\n",
    "\n",
    "    return Q, R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Q, R = gram_schmidt_qr(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.76781381 -0.59075334 -0.24793637]\n",
      " [-0.5817192  -0.48068519 -0.65615891]\n",
      " [-0.26844873 -0.64803723  0.71272929]]\n"
     ]
    }
   ],
   "source": [
    "print (Q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1.57007982e+00  -1.07526730e+00  -4.54661088e-01]\n",
      " [  2.87603125e-15   1.48523064e-01   1.10370436e-01]\n",
      " [ -5.90125834e-15   8.49331615e-15   2.89960049e-02]]\n"
     ]
    }
   ],
   "source": [
    "print (R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para verificarmos se a implementação está correta, devemos conferir se $QR = A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.20552897 -0.91334558 -0.42148593]\n",
      " [-0.91334558  0.5541108   0.19240567]\n",
      " [-0.42148593  0.19240567  0.07119534]]\n"
     ]
    }
   ],
   "source": [
    "print (Q.dot(R))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.20552897 -0.91334558 -0.42148593]\n",
      " [-0.91334558  0.5541108   0.19240567]\n",
      " [-0.42148593  0.19240567  0.07119534]]\n"
     ]
    }
   ],
   "source": [
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vou definir agora a minha implementação da $iteração QR$, algoritmo utilização para a obtenção dos autovalores da matriz $A$ (mais uma vez tomando como referência o livro do Justin Solomon):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def qr_iteration(A, steps=100):\n",
    "    for k in range(steps):\n",
    "        Q, R = gram_schmidt_qr(A)\n",
    "        A = R.dot(Q)\n",
    "    return np.diag(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.95951254 -0.15146022  0.0227828 ]\n"
     ]
    }
   ],
   "source": [
    "eigenvalues = qr_iteration(A)\n",
    "print (eigenvalues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usando o numpy como referência para conferir se a implementação está correta (a menos de um sinal talvez, como foi discutido em aula e apresentado no livro do Justin Solomon:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.95951254 -0.15146022  0.0227828 ]\n"
     ]
    }
   ],
   "source": [
    "w, _ = np.linalg.eig(A)\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}