Exemple algorithme QR pur pour le calcul des valeurs propres

algo_QR_pur.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:667dbdbb3dc845f9a3d6928eb522af241dfbbca538f793ba1f721665f17f9245"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import *"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def  householder(A): # computes a QR factorization upper triangular matrix R \n",
      "    \"\"\"matrices W, R pour la factorisation QR'\n",
      "    \n",
      "    Usage: \n",
      "    W, R = householder(A)\n",
      "\n",
      "    Argument:\n",
      "    A: matrice mxn avec m>=n\n",
      "    \n",
      "    Retour:\n",
      "    W: matrice de vecteurs de reflexions de householder\n",
      "    R: matrice R de la factorisation QR\n",
      "\n",
      "    Tire de: LN Trefethen & D Bau III, Numerical Linear Algebra, 1997 SIAM Philadelphia \n",
      "    \"\"\"\n",
      "    m,n = A.shape\n",
      "    ident = eye(m)\n",
      "    W = zeros((m,n))\n",
      "    R = copy(A)\n",
      " \n",
      "    for k in range(0,n):\n",
      "        e1=ident[k:m,k]\n",
      "        x = R[k:m,k]\n",
      "        W[k:m,k]= ( int(x[0]>=0) - int(x[0]<0) )*linalg.norm(x,2)*e1+x\n",
      "        W[k:m,k]=W[k:m,k]/linalg.norm(W[k:m,k],2)\n",
      "        R[k:m][:,k:n]=R[k:m][:,k:n]-dot(outer(2*W[k:m,k],W[k:m,k]),R[k:m][:,k:n])\n",
      " \n",
      "    return W, R\n",
      " \n",
      "def formQ(W): \n",
      "    \"\"\"Matrice Q de la factorisation QR\n",
      "    \n",
      "    Usage: \n",
      "    Q = formQ(W)\n",
      "    \n",
      "    Argument:\n",
      "    W: matrice mxn de reflexions obtenue de la d\u00e9composition Householder\n",
      "    \n",
      "    Retour:\n",
      "    Q: Matrice Q de la factorisation QR\n",
      "    \n",
      "    Note:\n",
      "    La matrice W se calcule avec W, R = householder(A)\n",
      "    \"\"\"\n",
      "    \n",
      "    m,n = W.shape\n",
      "    Q=zeros((m,m))\n",
      " \n",
      "    # compute the product QI = Q, column by column\n",
      "    for i in range(0,m):\n",
      "        Q[i,i]=1\n",
      "        for k in range(n-1,-1,-1): \n",
      "            Q[k:m,i]=Q[k:m,i]-2*dot(outer(W[k:m,k],W[k:m,k]),Q[k:m,i])\n",
      " \n",
      "    return Q"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 5\n",
      "A0 = random.random((n,n))-0.5\n",
      "A0 = dot(A0.T,A0) # on construit une matrice r\u00e9elle sym\u00e9trique\n",
      "print A0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 0.284666   -0.1626778  -0.07079496 -0.03470367  0.16145532]\n",
        " [-0.1626778   0.38794357  0.05822047  0.01157431 -0.12007222]\n",
        " [-0.07079496  0.05822047  0.3051367  -0.28920981 -0.01616102]\n",
        " [-0.03470367  0.01157431 -0.28920981  0.47465388 -0.09513347]\n",
        " [ 0.16145532 -0.12007222 -0.01616102 -0.09513347  0.11173422]]\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = A0\n",
      "# Algorithme QR 'pur'\n",
      "for k in range(0,20):\n",
      "    W, R = householder(A)\n",
      "    Q = formQ(W)\n",
      "    A = dot(R,Q)\n",
      "\n",
      "# affiche A avec 4 d\u00e9cimales\n",
      "for i in range(0,n):\n",
      "    for j in range(0,n):\n",
      "        print '{:+.4f}'.format(A[i,j]),\n",
      "    print '\\n'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "+0.6246 -0.0333 +0.0000 +0.0000 -0.0000 \n",
        "\n",
        "-0.0333 +0.6846 +0.0000 +0.0000 +0.0000 \n",
        "\n",
        "+0.0000 +0.0000 +0.1915 +0.0000 +0.0000 \n",
        "\n",
        "-0.0000 -0.0000 +0.0000 +0.0634 -0.0000 \n",
        "\n",
        "+0.0000 +0.0000 -0.0000 -0.0000 +0.0000 \n",
        "\n"
       ]
      }
     ],
     "prompt_number": 4
    }
   ],
   "metadata": {}
  }
 ]
}