nrla_naeherungsloesung_ueberbestimmtes_system.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Näherungslösungen von überbestimmten Gleichungssystemen\n",
    "\n",
    "- Gegeben: $A \\cdot x = b$ wobei $A \\in M(m \\times n), m > n$\n",
    "- Gesucht: Näherungslösung $\\bar{x}$ sodass $||b - A\\bar{x}|| \\leq || b - Ax||$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# input: equation system Ax = b\n",
    "A = matrix([\n",
    "        [1, 2],\n",
    "        [2, 3],\n",
    "        [3, 4]])\n",
    "b = vector([1, \n",
    "            2,\n",
    "            0])\n",
    "#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mit \"einfacher\" Pseudo-Inverse\n",
    "\n",
    "$$\\Sigma^{\\#} = (A^T \\cdot A)^{-1} \\cdot A^T$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}A{T} \\cdot A =  \\left(\\begin{array}{rr}\n",
       "14 & 20 \\\\\n",
       "20 & 29\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/plain": [
       "A{T} \\cdot A =  [14 20]\n",
       "[20 29]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}(A{T} \\cdot A)^{-1} =  \\left(\\begin{array}{rr}\n",
       "\\frac{29}{6} & -\\frac{10}{3} \\\\\n",
       "-\\frac{10}{3} & \\frac{7}{3}\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/plain": [
       "(A{T} \\cdot A)^{-1} =  [ 29/6 -10/3]\n",
       "[-10/3   7/3]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}A^{\\#} = (A{T} \\cdot A)^{-1} \\cdot A^{T} =  \\left(\\begin{array}{rrr}\n",
       "-\\frac{11}{6} & -\\frac{1}{3} & \\frac{7}{6} \\\\\n",
       "\\frac{4}{3} & \\frac{1}{3} & -\\frac{2}{3}\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/plain": [
       "A^{\\#} = (A{T} \\cdot A)^{-1} \\cdot A^{T} =  [-11/6  -1/3   7/6]\n",
       "[  4/3   1/3  -2/3]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Apinv = (A.T * A)^-1 * A.T\n",
    "show(LatexExpr(r\"A{T} \\cdot A = \"), A.T * A)\n",
    "show(LatexExpr(r\"(A{T} \\cdot A)^{-1} = \"), (A.T * A)^-1)\n",
    "show(LatexExpr(r\"A^{\\#} = (A{T} \\cdot A)^{-1} \\cdot A^{T} = \"), Apinv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\bar{x} = \\left(-\\frac{5}{2},\\,2\\right)</script></html>"
      ],
      "text/plain": [
       "\\bar{x} = (-5/2, 2)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = Apinv * b\n",
    "show(LatexExpr(r\"\\bar{x} =\"), x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}A \\cdot \\bar{x} = \\left(\\frac{3}{2},\\,1,\\,\\frac{1}{2}\\right)</script></html>"
      ],
      "text/plain": [
       "A \\cdot \\bar{x} = (3/2, 1, 1/2)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}b =  \\left(1,\\,2,\\,0\\right)</script></html>"
      ],
      "text/plain": [
       "b =  (1, 2, 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(LatexExpr(r\"A \\cdot \\bar{x} =\"), A*x)\n",
    "show(LatexExpr(r\"b = \"), b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mit Pseudoinverse ermittelt mit SVD\n",
    "\n",
    "Singulärwertzerlegung $A = U \\cdot \\Sigma V^{T}$\n",
    "\n",
    "Die Pseudoinverse ist definiert als:\n",
    "$$ A^{\\#} = V \\cdot \\Sigma^{\\#} \\cdot U^{T} $$\n",
    "mit $\\Sigma^{\\#} = \\begin{pmatrix}D^{-1} & 0 \\\\ 0 & 0 \\\\ \\end{pmatrix} \\in M(n \\times m)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
       "-0.338098165838 & 0.847952217752 & 0.408248290464 \\\\\n",
       "-0.550649318286 & 0.173547289246 & -0.816496580928 \\\\\n",
       "-0.763200470734 & -0.50085763926 & 0.408248290464\n",
       "\\end{array}\\right) \\left(\\begin{array}{rr}\n",
       "6.54675563644 & 0.0 \\\\\n",
       "0.0 & 0.37415322624 \\\\\n",
       "0.0 & 0.0\n",
       "\\end{array}\\right) \\left(\\begin{array}{rr}\n",
       "-0.569594837763 & -0.821925617556 \\\\\n",
       "-0.821925617556 & 0.569594837763\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/plain": [
       "[ -0.338098165838458  0.8479522177516872  0.4082482904638629]\n",
       "[-0.5506493182859842 0.17354728924565255 -0.8164965809277263]\n",
       "[-0.7632004707335103 -0.5008576392603808  0.4082482904638632] [ 6.546755636442668                0.0]\n",
       "[               0.0 0.3741532262404972]\n",
       "[               0.0                0.0] [-0.5695948377626012 -0.8219256175556253]\n",
       "[-0.8219256175556253  0.5695948377626012]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "U, S, V  = matrix(RDF, A).SVD()\n",
    "show(U,S,V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
       "-1.83333333333 & -0.333333333333 & 1.16666666667 \\\\\n",
       "1.33333333333 & 0.333333333333 & -0.666666666667\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/plain": [
       "[-1.8333333333333306 -0.3333333333333316  1.1666666666666647]\n",
       "[ 1.3333333333333308 0.33333333333333204  -0.666666666666665]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Si = matrix(RDF, A.ncols(), A.nrows())\n",
    "for i in range(S.ncols()):\n",
    "    Si[i,i] = S[i,i]^-1 if S[i,i] != 0 else 0.0\n",
    "Apinv = V *  Si * U.T\n",
    "show(Apinv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\bar{x} = \\left(-2.5,\\,2.0\\right)</script></html>"
      ],
      "text/plain": [
       "\\bar{x} = (-2.499999999999994, 1.999999999999995)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = Apinv * b\n",
    "show(LatexExpr(r\"\\bar{x} =\"), x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}A \\cdot \\bar{x} = \\left(1.5,\\,1.0,\\,0.5\\right)</script></html>"
      ],
      "text/plain": [
       "A \\cdot \\bar{x} = (1.499999999999996, 0.9999999999999973, 0.4999999999999982)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}b =  \\left(1,\\,2,\\,0\\right)</script></html>"
      ],
      "text/plain": [
       "b =  (1, 2, 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(LatexExpr(r\"A \\cdot \\bar{x} =\"),A*x)\n",
    "show(LatexExpr(r\"b = \"), b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mit QR Zerlegung\n",
    "\n",
    "QR Zerlegung $A = QR$\n",
    "\n",
    "Näherungslösung wird berechnet durch lösen von:\n",
    "$$R \\bar{x} = Q^T b$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}Q = \\left(\\begin{array}{rrr}\n",
       "-0.267261241912 & 0.872871560944 & 0.408248290464 \\\\\n",
       "-0.534522483825 & 0.218217890236 & -0.816496580928 \\\\\n",
       "-0.801783725737 & -0.436435780472 & 0.408248290464\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/plain": [
       "Q = [-0.2672612419124243  0.8728715609439699  0.4082482904638629]\n",
       "[-0.5345224838248488 0.21821789023599203 -0.8164965809277263]\n",
       "[-0.8017837257372732 -0.4364357804719846  0.4082482904638632]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}R = \\left(\\begin{array}{rr}\n",
       "-3.74165738677 & -5.34522483825 \\\\\n",
       "0.0 & 0.654653670708 \\\\\n",
       "0.0 & 0.0\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/plain": [
       "R = [-3.7416573867739413  -5.345224838248488]\n",
       "[                0.0  0.6546536707079786]\n",
       "[                0.0                 0.0]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "try:\n",
    "    # try rational field first\n",
    "    Q,R = A.QR()\n",
    "except TypeError:\n",
    "    # if it doesn't work, use real double field\n",
    "    Q,R = matrix(RDF, A).QR()\n",
    "    \n",
    "show(LatexExpr(r\"Q =\"), Q)\n",
    "show(LatexExpr(r\"R =\"), R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}R x =  \\left(-3.74165738677 \\, x_{0} - 5.34522483825 \\, x_{1},\\,0.654653670708 \\, x_{1},\\,0.0\\right)</script></html>"
      ],
      "text/plain": [
       "R x =  (-3.7416573867739413*x0 - 5.345224838248488*x1, 0.6546536707079786*x1, 0.0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}Q^T b =  \\left(-1.33630620956,\\,1.30930734142,\\,-1.22474487139\\right)</script></html>"
      ],
      "text/plain": [
       "Q^T b =  (-1.3363062095621219, 1.309307341415954, -1.2247448713915896)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Xsym = vector([var(\"x{}\".format(i)) for i in range(A.ncols())])\n",
    "RX = R * Xsym\n",
    "QTb = Q.T * b\n",
    "show(LatexExpr(r\"R x = \"), RX)\n",
    "show(LatexExpr(r\"Q^T b = \"),QTb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No solution found\n"
     ]
    }
   ],
   "source": [
    "sol = solve([RX[i] == QTb[i] for i in range(A.nrows())], *Xsym)\n",
    "if sol:\n",
    "    for s in sol:\n",
    "        show(s)\n",
    "else:\n",
    "    print \"No solution found\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 7.0",
   "language": "",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}