{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nonnegative matrix factorization\n",
    "\n",
    "A derivative work by Judson Wilson, 6/2/2014.    \n",
    "Adapted from the CVX example of the same name, by Argyris Zymnis, Joelle Skaf and Stephen Boyd\n",
    "\n",
    "## Introduction\n",
    "\n",
    "We are given a matrix $A \\in \\mathbf{\\mbox{R}}^{m \\times n}$ and are interested in solving the problem:\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{minimize}   & \\| A - YX \\|_F \\\\\n",
    "    \\mbox{subject to} & Y \\succeq 0 \\\\\n",
    "                      & X \\succeq 0,\n",
    "    \\end{array}\n",
    "where $Y \\in \\mathbf{\\mbox{R}}^{m \\times k}$ and $X \\in \\mathbf{\\mbox{R}}^{k \\times n}$.\n",
    "\n",
    "This example generates a random matrix $A$ and obtains an\n",
    "*approximate* solution to the above problem by first generating\n",
    "a random initial guess for $Y$ and then alternatively minimizing\n",
    "over $X$ and $Y$ for a fixed number of iterations.\n",
    "\n",
    "## Generate problem data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "\n",
    "# Ensure repeatably random problem data.\n",
    "np.random.seed(0)\n",
    "\n",
    "# Generate random data matrix A.\n",
    "m = 10\n",
    "n = 10\n",
    "k = 5\n",
    "A = np.random.rand(m, k).dot(np.random.rand(k, n))\n",
    "\n",
    "# Initialize Y randomly.\n",
    "Y_init = np.random.rand(m, k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perform alternating minimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 1, residual norm 2.766300564135502\n",
      "Iteration 2, residual norm 0.5840356930600721\n",
      "Iteration 3, residual norm 0.3356679970549085\n",
      "Iteration 4, residual norm 0.18670276027770083\n",
      "Iteration 5, residual norm 0.12819921698143966\n",
      "Iteration 6, residual norm 0.09295501592922492\n",
      "Iteration 7, residual norm 0.06766021043574907\n",
      "Iteration 8, residual norm 0.04958204907945361\n",
      "Iteration 9, residual norm 0.03897402158866238\n",
      "Iteration 10, residual norm 0.02979328283505179\n",
      "Iteration 11, residual norm 0.022938564327729952\n",
      "Iteration 12, residual norm 0.021943924920767337\n",
      "Iteration 13, residual norm 0.01810297853945281\n",
      "Iteration 14, residual norm 0.014551161988556204\n",
      "Iteration 15, residual norm 0.014039687334395924\n",
      "Iteration 16, residual norm 0.009354606824469416\n",
      "Iteration 17, residual norm 0.008643141637584189\n",
      "Iteration 18, residual norm 0.007278100007476402\n",
      "Iteration 19, residual norm 0.008486679700021057\n",
      "Iteration 20, residual norm 0.008827511916396866\n",
      "Iteration 21, residual norm 0.008396764193205366\n",
      "Iteration 22, residual norm 0.005265185332845983\n",
      "Iteration 23, residual norm 0.006931929503816392\n",
      "Iteration 24, residual norm 0.007356156596477946\n",
      "Iteration 25, residual norm 0.0039053948996930054\n",
      "Iteration 26, residual norm 0.003989885269615319\n",
      "Iteration 27, residual norm 0.002920361405226024\n",
      "Iteration 28, residual norm 0.007779246694466739\n",
      "Iteration 29, residual norm 0.007339011292898449\n",
      "Iteration 30, residual norm 0.005008539285258121\n"
     ]
    }
   ],
   "source": [
    "# Ensure same initial random Y, rather than generate new one\n",
    "# when executing this cell.\n",
    "Y = Y_init \n",
    "\n",
    "# Perform alternating minimization.\n",
    "MAX_ITERS = 30\n",
    "residual = np.zeros(MAX_ITERS)\n",
    "for iter_num in range(1, 1+MAX_ITERS):\n",
    "    # At the beginning of an iteration, X and Y are NumPy\n",
    "    # array types, NOT CVXPY variables.\n",
    "\n",
    "    # For odd iterations, treat Y constant, optimize over X.\n",
    "    if iter_num % 2 == 1:\n",
    "        X = cp.Variable(shape=(k, n))\n",
    "        constraint = [X >= 0]\n",
    "    # For even iterations, treat X constant, optimize over Y.\n",
    "    else:\n",
    "        Y = cp.Variable(shape=(m, k))\n",
    "        constraint = [Y >= 0]\n",
    "    \n",
    "    # Solve the problem.\n",
    "    # increase max iters otherwise, a few iterations are \"OPTIMAL_INACCURATE\"\n",
    "    # (eg a few of the entries in X or Y are negative beyond standard tolerances)\n",
    "    obj = cp.Minimize(cp.norm(A - Y*X, 'fro'))\n",
    "    prob = cp.Problem(obj, constraint)\n",
    "    prob.solve(solver=cp.SCS, max_iters=10000)\n",
    "\n",
    "    if prob.status != cp.OPTIMAL:\n",
    "        raise Exception(\"Solver did not converge!\")\n",
    "    \n",
    "    print('Iteration {}, residual norm {}'.format(iter_num, prob.value))\n",
    "    residual[iter_num-1] = prob.value\n",
    "\n",
    "    # Convert variable to NumPy array constant for next iteration.\n",
    "    if iter_num % 2 == 1:\n",
    "        X = X.value\n",
    "    else:\n",
    "        Y = Y.value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Output results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEYCAYAAABcGYHrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAHT1JREFUeJzt3WtwXOd93/Hf/+wuQIAiAIKkZVUyTUKMZNmyY/Mity/8QhZoN3WctA5l2e30nUTa4xeddGRSyrR2Om1HotROW2emKcmmM51OkuriTJMX6dikNBnHaUfmxWlatbqY0MWSo4sNgKREkAB2/31xngUOFrvA4mB3zy7O9zODwe7BXv7cI53fPuc5z/OYuwsAgCjrAgAA3YFAAABIIhAAAAGBAACQRCAAAAICAQAgiUAAAAQEAgBAEoEAAAiKWRewFtu3b/ddu3ZlXQYA9JTz58//3N13rPa4ngqEXbt26dy5c1mXAQA9xcxea+ZxnDICAEgiEAAAAYEAAJBEIAAAAgIBACCJQAAABAQCAEBSTgLhz19+V//2zEtZlwEAXS0XgfDcxKS+88zLqlRYPxoAGslFIAwPlFRx6b3Z+axLAYCulZtAkKRLV+cyrgQAulcuAmGoGggzBAIANJKLQBgmEABgVbkIhJFBAgEAVpOLQKCFAACrIxAAAJJyEgiDfQUVI9M0VxkBQEO5CAQz0/BAiRYCAKwgF4EgScODJV0mEACgofwEAi0EAFhRrgJhemY26zIAoGvlKhBoIQBAY7kJhJGBEnMZAcAKchMIwwMlXbk+zxTYANBAbgJhaKAkd+nKNabABoB6chMI1dHKdCwDQH25CwQ6lgGgvtwEwshgnyQCAQAayU0g0EIAgJURCAAASTkMBGY8BYD6chMIm0qR+goRE9wBQAO5CQQz0/Ag01cAQCO5CQSJ+YwAYCUEAgBAUg4DgU5lAKgvd4FACwEA6stdIHCVEQDUl7tAuHJ9XmWmwAaAZXIXCJJoJQBAHbkMhGkCAQCWyVUgjAwynxEANJKrQGCCOwBojEAAAEjKayBcZRlNAKiVq0AYooUAAA3lKhA2lQraVIoIBACoI1eBIDF9BQA0QiAAACTlNBCY8RQAlstlINBCAIDlchgIfcxlBAB15DAQaCEAQD25DIT3Z8uaK1eyLgUAukoOA6EoicFpAFArd4EwMtgniUAAgFotDQQz+2wrX68dmOAOAOorpn2imX1J0mjN5iOSDqyrojZjPiMAqC9VIJjZf5A0Jmla0mTYPCpppEV1tQ3LaAJAfWlbCOfd/Wu1G83sgWaebGYjkg6HuwcknXb3kylrWZOFZTQZrQwAS6QNhMkG2083+fyH3f1Y9Y6ZXTQzdSIU6EMAgPrSdipfMLPPmtkuMxuq/kg6ttoTQ+tgrGbziWae2wp9xUiDfQUCAQBqpG0hjCs+iFe5JAu/v77Kc0cljZvZmLtPhG3TWh4SbcNoZQBYLm0gjEja6u6XkhvN7NHVnhhCYGvN5oOSztR7vJkdVuhv2LlzZ6piaxEIALBc6lNGtWEQPLLWFwqnkMbV4JSRu5909/3uvn/Hjh1rffm6hgZKukSnMgAskTYQ3Mx21dne1FVGNU5JusfdL6SsZc1GaCEAwDJpTxl9TdKnzEySqv0A2yTtlvSvmn0RMzsq6UQnw0DilBEA1JM2EMYkPaS4M7jKJB1t9gXM7JDiU09nwv3x6u12IxAAYLm0gXDM3Z+p3Whmv2jmyWY2rvhqozNmVr26qGHHcqsND5Q0M1fW9fmy+ouFTrwlAHS9tH0Iw2b2ydqN7v7j1Z4YOpFPK75sdUrSxfDTsctORwYZnAYAtdIGwlfqbQyD01bk7tPubnV+7k1Zy5oNMZ8RACyTNhCe0GJnctLhOtu6DtNXAMByafsQDkp61MwmtNixbJLu0RquMsoKgQAAy6UNhP2SHtPySe66fvpriRlPAaCeVl9lNF3vwd2GZTQBYLlUfQjVMAiznH6y2plcLyS60dCmOAcJBABYlHpN5bBq2rSkZyVNmdkTLauqzYqFSDf0FwkEAEhIFQhm9k3Fq5xF7j7q7gVJT5rZg60tr30YrQwAS6VtIUy4+3eTG8L9ejOgdiVmPAWApVLPdtpge1NTV3QDZjwFgKXSBsKttaOSw3TYd623oE7hlBEALJX2stOTkp41M1c8FmFU8RiEfa0qrN0IBABYKlUghNXS9pvZbyielG5Zn0K3Gx4kEAAgKW0LQdJCR3JPGh4o6fp8RdfmytpUYgpsAGiqD8HM1rxWcrdjPiMAWKrZFsKRsPjNhOJJ7KpXGVVvV/sQ5O5dP7mdtDQQbhzalHE1AJC9ZgPhTKMDfehHeEzxIjcdW9NgvWghAMBSzV52uuyUUZjH6ElJT0p6xN0PuPurrSyunRYCgcFpACCpyUCoXRrTzL4k6VVJw5L2uPvjrS+tvarLaE7TQgAASWu8yigMRnta8XiDY+5+qi1VdQCnjABgqaZHKpvZA4pnN3VJY70cBpK0ZROBAABJzV52+n1Jj0o65O6fDwPT6j2uZy5PLUSmLZuKukwgAICk5k8Z7VfcsbzNzO5XfLmptHSSuz2SHpD0cOvKay+mrwCARc0GwslmOo7NbGyd9XTUyGBJ01dnsy4DALpCs30IJ5p83LG0hWSBFgIALGr2stNXWvm4bkEgAMCi1GsqbwRxIMxnXQYAdIWcB0KfLs3Myr3RAnAAkB85D4SS5squmbly1qUAQOZyHwgSg9MAQCIQJBEIACC1OBDM7MFWvl67MeMpACxqamCamX2vmYcpnvSuJxbIkZjxFACSmh2pvE2rDzozSUfXV05nccoIABY1GwgP1K6JUE9YZrNnDIVAYII7AEi5QE4tMxsOS2n21AX9W/qLMqOFAADSOjuVwzKauyRtlXRe0n0tqKljosg0tInpKwBAWuOKaVVmdo+kp7TYIrBw+4EW1dUx8YynBAIApG0hjLv7qLtvk3TY3UcljSleUa2nMMEdAMTSBsK5xO2tktRoFbVuRyAAQCx1H4KZfSnc3Gpmvxxu711/SZ01NFDiKiMAUPpAmJD0W6FD+aSk74ZLTm9tUV0dM0ILAQAkpexUDpeh7k9s2mNmn2pmrEK3GR4oaXpmTu4uM1v9CQCwQbVsLiN3/3FoMfSU4YGSyhXX+7NMgQ0g31IFQhh/sOxH0vEW19d2TF8BALFUp4wUX17qiscfJE2tr5zOS854evPIQMbVAEB20p4yetrdC+4eVX8Uz3Q63sLaOmJ4kBYCAEjpA2HZiOTQobx7feV03uIpo9mMKwGAbKUKhBUGoY2uo5ZM0IcAALG0cxnVWzBnTNLT6yun8wgEAIil7VTeJukRLZ27aMLdX1l/SZ11Q39RhcgIBAC5lzYQjrn7My2tJCNmxnxGAKD0fQh1wyCMReg5wwNMgQ0ALRupHJxq8et1xBAtBABo7pSRmVXUY8tjrsXwQEmXrnLZKYB8a7aFUB2IVnD3gqTPSdpTs+1A2N5zmPEUAJoPhNqBaMO1VxS5+wX1aCuiOuMpAORZU4FQZyBaowFoI+srJxvDYZGcSqUn8wwAWiJtp/IeM/tkckO4f2D9JXXe8EBJFZfem53PuhQAyEzaBXIeMrPvm9luxYPTxhSvonZPK4vrlOSMp0ObShlXAwDZSDswTe7+OTPbq3iW04m1DFQzszFJhyRtc/djaWtoleSMpx/KuBYAyErqQJAWOpIvVO+b2f3u/h9Xeo6ZjSvua+ia9ZeZzwgAmh+H8LuSnnL3Z8P9epPbmeLWwoqB4O5nwmscUJd0QhMIANB8C6F2ZbRtkmpP9Ziko+uuKAMEAgA0GQju/rWaTQ+EBXGWMLNftKSqDiMQACD9ZadT1ctOzWzIzB40swfrhcR6mdlhMztnZufefffdVr+8JGmwr6BSgSmwAeRb2kB4SPGlppL0jOJTSM+Y2YMtqSrB3U+6+353379jx45Wv7ykxSmwmfEUQJ6lvcrotLv/URiHsM/dD0gLl5P2pKEwWhkA8ir1KaPwe1xLl83s2bkfWCQHQN6lbSHsM7Otiq80OixJZnaPGs9x1PVGBkr6+XtMgQ0gv9KumPa44oP/EXd/NoTB3maea2Z7zey44pHK42Z2PIx4zhQtBAB5t56Ryk9I+rKZ7XL3Z8xsspmrjBKjmzOfsiIp7lSmhQAgv1K1EEKL4FnFI5OrHcnTZvbZVhXWacMDJV25Ps8U2AByK22n8sFwKejXFUYxhwVztrassg4bGijJXbpyjSmwAeRT2kD4UYPtPfv1emSwTxKjlQHkV9pAuMvMtoTbLklmtkvSXS2oKRNMXwEg79J2Kj8i6cdmNiVJZjaieKGcnlwgR1oMhOkZOpYB5FPaFdMuKV5G8ze0uFradCsL6zRaCADybr0L5Hw3ed/M7tcq6yF0KwIBQN413YcQZjW938zubvR3ddEqaGs1MkggAMi3pgLBzIYVDyY7KemMmf29sP1BMzsb1kGYUpesgJbGplJBfcVIl5jxFEBONXvK6GFJx939VJjR9KiZbZP0FUnnFE+Bfbb2FFKvYfoKAHnWdB+Cu58KvyfM7ISkh9x9f9sqywCBACDPmu1DWLI0Zpiz6MnWl5OtEQIBQI41Gwj1RiAv29aOFdM6iRYCgDxr9pTRfWYmLR1rsM/MkusfjKiHrzKS4kB44a0rWZcBAJloNhD2qf7EdQdr7p9fXznZYhlNAHnWbCA85u4PrfYgM3t0nfVkqjoFdrniKkSWdTkA0FHN9iGcaPHjulJ1cBqtBAB51FQghLUOWva4bsX0FQDyLO301xvS4oynBAKA/CEQEmghAMgzAiGBQACQZwRCwjAzngLIMQIhYetgnwb7Cvrhy+9mXQoAdByBkFAqRPrG3Xv0veff1p8TCgByhkCocf9nduvD2wb123/yvGbnK1mXAwAdQyDU6C8W9K1f/aguvvu+/vP/eDXrcgCgYwiEOu6540bdffsO/btnXtY7l69lXQ4AdASB0MC3vvgxzc5X9Oh/fyHrUgCgIwiEBnZv36z7P7Nbf/TjN3X+tcmsywGAtiMQVvCNu/fog0Ob9K0/fl7lSr01ggBg4yAQVrC5v6jf+sIdev5nl/Vfz76edTkA0FYEwiq++Imb9Ondo3r8ey9q6v3ZrMsBgLYhEFZhZvrtX/uYrlyb178+/WLW5QBA2xAITbjjpiH9w7/5Yf3Bc6/r+Z9dyrocAGgLAqFJvzl+m0YG+/TtP35e7nQwA9h4CIQmDQ+WdPTzt+vca1P6b3/5ZtblAEDLEQhr8OX9H9Iv3zKsR/70Bb13fT7rcgCgpQiENYiiuIP5nSvX9TvPvJx1OQDQUgTCGn1q51bdu+8W/d4PX9FP3nkv63IAoGUIhBSO/u2PaKBU0Lf/5P9ovswU2QA2BgIhhR1b+vXw37lDf/GTX+ibT/8V01oA2BCKWRfQq/7+p3dq6uqsHv/ei4rM9PihTyiKLOuyACA1AmEdvnH3Hs2XXf/mzEsqRqZHvvRxQgFAzyIQ1ukfjf+SypWKvvPsTxRFpn/5d+8kFAD0JAKhBX7z4G2ar7j+/Z9dVCGS/vmv3ykzQgFAbyEQWsDM9M3P365yxXXiBxMqRpG+/cWPEgoAegqB0CJmpod+5SOar7h+74evKDLTP/3VOwgFAD2DQGghM9M/+cIdKldc/+kvXlGxYHr4Vz5CKADoCQRCi5mZvv3Fj6pccZ38wYQKkeno528nFAB0PQKhDcxM/+zXPqayu373zy6qFJn+8eduz7osAFgRgdAmUWT6F79+p8pl13ee/YkuvD6tf/DpnRr/6I0qFRggDqD7EAhtFIXBaru2b9Z/+Z+v6uu/f0Hbb+jXvftv0VcP7NTObYNZlwgAC6yXVv/av3+/nzt3LusyUilXXD946V39wY9e17MvvKNyxfWZX9qur961U+N33Ki+Iq0GAO1hZufdff+qjyMQOu+tS9f05Lmf6omzP9Wb0zPafkOfDu37kL5y4EPatX1z1uUB2GAIhB5Qrrh+8PK7+sPnXtczodVwYNdW/a2xbbpr9zZ9aueINvdzVg/A+hAIPebty9f05Nmf6vv/9209/7NLqrhUiEx33jysu3Zt1V27t+nArq0aGezLulQAPYZA6GFXrs3pwuvTOvvKpH70yqT+8o1pzc7HC/HcfuMWHdi9VQd2jerjNw/rw9s2q8BkegBWQCBsINfmyvqrNy7p7KuTeu6VSV14bUrvXZ+XJPUXI+35wA26/cYtuu2DW3T7B7fo9hu36KbhTQyGAyCJQNjQ5ssVvfDWFf2/v76sl96+ohfffk8vvXVFb12+tvCYLf1F3fbBLbrtxi3a84EbdMvWgfhnZFBDA0XCAsiRZgOBHsseVCxEuvPmYd158/CS7Zeuzumld67oxbfCz9tX9Kf/+691aWZuyeO29Bd1cwiIm0cGdMvWQd0cbm/f0q+RgZIG+wqEBpAzmQWCmR2VNCFpVJLc/WRWtWwUw4MlHdg1qgO7Rhe2ubumrs7pzakZvTF1VW9Oz+iNqerPVT03Makr4fRTUl8h0shgKfz0aWSgpK2DfYv3B0sa2lTS0EAx/C5paFNRWzaVGFMB9KhMAsHMjks66+5PV++b2aHqfbSOmWl0c59GN/fp47cM133MpZm5OCymZjR1dVZTV+c0fXVO01dnNX11TlNXZ/X65FX9rzemNXV1bqGDu5GBUmFJUGzuL2pzX0GDfUVt7i9ooK+gzX1FDfYVtLk//O4ralOpoMgkmRSZKTKTmRRZ/O+IzGSK/1YqmkqFSH2FSP3FSH3FKL5fjFSMjNYNkEJWLYTD7n4scf+0pGOSCIQMDA+UNDwwrI/9jfqBkeTuujZX0fTMrK5cm9flmTldvjanyzPz4fecLie2X5qJt711aUbvXy/r6uy83p8trxoq62EWt3D6itHC71IhUqlgC6HRVwjbipH6wvbITAo5YopDyMLrVe9L8eXA/cVIm0oF9Rcj9RcL6i9FC7c3lcK2YqQokkxxsCVfrxpusvjvkcWvGwdffLsaiNXbUfibh/1QrkgVd1Xc5V69rXA/vl3tInR3VXsL3aWFe8kuxIV/uynxUcR1Jz6DvkKkYvjMqp9psVDdHrZFkSruKofaypX4dqXiidta2FatKy7JQ41x3ckyS1G08GWg+oWgVLCFz67T5ssVXZ+v6NpceeH3bLki93g/S/U+z7A1/HdQjOL6F39HKhSWbu/Uv63jgWBme+tsnpQ03ulasHZmpoG+ggb6BnTT6vnR0Fy5oquzISBCUFybqyweyLT0IOc1B735ckWz4X/GuXJFs/PxT/X29cS2+bJrrhxvm6s+plzR3LxrZmZu4TmVxEFTiQNS/HvxIFouu67PV8JPWXPl3rkwY6Na+BIQQj6qewBttJ/ioK0GrpkpihKt1PD67lpy4L8+X9F8pTP7PjLpd766V1/4xE1tfZ8sWgijigMgabrRg83ssKTDkrRz5842loVOKhUiDQ9EGh4oZV3KupUrruvzZV2fWwyJ6kFjIcy0GGrVgKks3K/5hl/xhdvlSvxtv5z4m9niwaoQLZ5Oqx7ULPG7+u1U0kKLJP61dPuy1kOiAbEkDCu+ELBzFdfcfEXzlYpmy4u358Lf4/oWa6veLkSmKDIVEvUnW2RSaJXV1C9p8b3L4T1DyM8mQn6uHB+o62VCo5jw0GJJfgFZuK14u0naVFraAqy2EpO/+4rRwvvUfn61LZ9KaOmVK3HN1c+3XPFwf3H7ng/c0NR/j+uRRSCM1Nk2KUlmNuLuS8IhdDaflOLLTttfHrA2hcg02FcUg8jR67K4HGRa4cqihNr7AIAOyyIQJrW8lTAiSbWtAwBA53Q8ENz9gpb3GYxKOtPpWgAAi7IaQfSkmR1K3D8o6URGtQAAlFEguPsRSWNmdiiMWL7IoDQAyFZmU1e4+2NZvTcAYDkmnQEASCIQAABBT62HYGbvSnot5dO3S/p5C8tB67GPuh/7qPvV20cfdvcdqz2xpwJhPczsXDMLRCA77KPuxz7qfuvZR5wyAgBIIhAAAEGeAoEV2bof+6j7sY+6X+p9lJs+BADAyvLUQgAArIBAAABIIhAAAEFmcxl1Spg8b0JhEZ6wAhsyYmZjkg5J2ubux+r8nf2VMTMbUVi2VtIBSadr9wP7KVthH3053L1Vkmr/f0qzjzZ0IJjZcUlnqzOpmtlxMzvEzKrZMLNxxYsh3drg7+yv7vBw8uBiZhfNbOGAwn7qCsclHasuKmZm583saHXS0LT7aKOfMjpc8wGclnQkq2Lyzt3PhP3RaGU89lfGwjfPsZrNJyQlv32yn7K3X9J44v6E4tZcVap9tGEDwcz21tk8qaUfIroE+6trjEoaD6f2qqYVQoL91B3cfV/NAX+v4oP+uvbRRj5lNKr4Q0hizebuxf7qAu4+IWlrzeaDWlzilv3UZUJfwZlEH0HqfbRhWwiKz1XXmpQWmsXoLuyvLhQ++3EtnjJiP3UJMxsxs2rn/8XEn1Lvo43cQphW6F1PqL2P7sH+6k6nJN3j7hfCffZTlwgdytWO/qfM7Cl3v1fr2EcbuYUwqeVJOSItfJDoLuyvLhNORZxIhIHEfspcaBkcrdl8WvHl3NI69tGGDYTwH3HtP35Ui+dC0UXYX93FzA5JuuDuZ8L9cYn91CX2Szpec/pn4fZ69tGGDYTgyfAfdtVBxZfQoTuxv7pAOPiPSjpnZmPhiqODiYewnzIUQvpYzbf9g5IeS9xPtY82/GynidF6Y5KmGVGZnXA53H1abNo+LemJ5CkJ9le2wrfOqTp/ejqcn64+jv2UocSIf0naJukX1UFpiceseR9t+EAAADRno58yAgA0iUAAAEgiEAAAAYEAAJBEIAAAAgIBACCJQAAABAQCAEASgQAACAgEtI2Z7Q1ruXp1zdesa0oysxNm1tY5eMJcQMfNbMrMTtf8bW+YtngqrIHbjvcfD5/9U+14fWwsBALaxt0vhMXaJxTPWVQ718rh+s9svQbv9VT4aRt3nwifwSOS9ifrCJ/PvZIeSS5q3+L3PyMmnkOTCARkaV+W7+XuZ6rTO3fAtKQHJJ2os2rVRJvfu3Y5RaAuAgGZCKePOrLSViffayVhUfQzanOrBEhrIy+hiS4V5ts/IGlvOFgvTM0bvj0/LOlseMxpdz8TnnNc0jnFq0Pdp/g01NNhWu2x8PILz1npvcJzTkmadPeDifc+rMVv7GPV01yJ95/Q4imYg5IurnHq53slvWJmh0JA1PtsjlfrCjWdUjzV8VZ3n25Qy17FrZAz4faopH3ufiT8fSQxP/6oJNXWneazX8O/G73A3fnhp60/ihcAP1qz7ZCkpxo8diRx/3z1vuKD9UXFq0PtlbQ38ZhDiedMNfle44oPesveK/G8EzX3LyoOCoU6vMnP4HDytuI1B6r/rkN16j1ds83r1LasFknjicdUl1U8pDi4kq93IllT2s+en431wykjdI3wDXbal64EdU7Sl8PtSUkT7j7tcYdsdWGde73m22qd8/TNvLeS7x1e83DitabD9onax66Fx9/Mzyn+5l9PM6/bqJaJmsdUT5Ul10WW4tNWC1c2reOzxwbCKSN0kzFpcf3e4CktP8jVmqyeDtJiB+pog8eu9N71Ol+nw9+qB8BlHcBmNpIiHI5IulizzOFa1euMnmxwu95zk6GZ9rPHBkIgIHNhOUApHKR87Vf+nFfcSrgQXu9U+L3sQF19r+o364QJ1e94HlEbrgJy9wkzO6a4ldCWS05XUfvvSvvZYwPhlBGyUl3rVYrPg0+EUzTTiYCoDuza2+hFwt9GE2GQ/NZb/ba77L1qXye890jNex9SvJZwMlTSXq10a533fEzLv6nX1qsV/v2r1ZL8+1jN344occoozWePjYdAQNtURyorPhjdlxypHA7gE2GgVvJgdY+kI2Z2KByQx9z9QjiVcUTSuJkdrR74w+s8GbaNS9qv+Hr/IwrfgOu9VzjQHdPSwWL7Eu99OLz3vTWPH6u+f2J08fHkgbTmMxgLo4QPNxiNfG/thhBYT5vZ4cQppWlJpxIH6RVrCfXvDZ/DmKQHqp9p2A+nffnVUWv67LHxmMdXEAAAco4WAgBAEoEAAAgIBACAJAIBABAQCAAASQQCACAgEAAAkggEAEBAIAAAJBEIAIDg/wPpi0sDT5HS4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original matrix:\n",
      "[[1.323426   1.11061189 1.69137835 1.20020115 1.13216889 0.5980743\n",
      "  1.64965406 0.340611   1.69871738 0.78278448]\n",
      " [1.73721109 1.40464204 1.90898877 1.60774132 1.53717253 0.62647405\n",
      "  1.76242265 0.41151492 1.8048194  1.20313124]\n",
      " [1.4071438  1.10269406 1.75323063 1.18928983 1.23428169 0.60364688\n",
      "  1.63792853 0.40855006 1.57257432 1.17227344]\n",
      " [1.3905141  1.33367163 1.07723947 1.67735654 1.33039096 0.42003169\n",
      "  1.22641711 0.21470465 1.47350799 0.84931787]\n",
      " [1.42153652 1.13598552 2.00816457 1.11463462 1.17914429 0.69942578\n",
      "  1.90353699 0.45664487 1.81023916 1.09668578]\n",
      " [1.60813803 1.23214532 1.73741086 1.3148874  1.27589039 0.40755835\n",
      "  1.31904948 0.3469129  1.34256526 0.76924618]\n",
      " [0.90607895 0.6632877  1.25412229 0.81696721 0.87218892 0.50032884\n",
      "  1.245879   0.25079329 1.25017792 0.72155621]\n",
      " [1.5691922  1.47359672 1.76518996 1.66268312 1.43746574 0.72486628\n",
      "  1.97409333 0.39239642 2.09234807 1.16325748]\n",
      " [1.18723548 1.00282008 1.41532595 1.03836298 0.90382914 0.38460446\n",
      "  1.213473   0.23641422 1.32784402 0.27179726]\n",
      " [0.75789915 0.75119989 0.99502166 0.65444815 0.56073096 0.341146\n",
      "  1.02555143 0.24273668 1.01035919 0.49427978]]\n",
      "Left factor Y:\n",
      "[[ 7.56475742e-01  3.42102372e-01  8.40426641e-01  7.02845111e-01\n",
      "   4.38002833e-03]\n",
      " [ 6.36189366e-01  8.27831861e-01  5.28165827e-01  5.60609403e-01\n",
      "   3.34595403e-02]\n",
      " [ 5.54834858e-01  6.37954560e-01  8.01726231e-01  1.96879041e-01\n",
      "   3.74736667e-02]\n",
      " [ 2.72955779e-01  9.53749151e-01  6.14934798e-02  9.81276972e-01\n",
      "  -4.26647247e-05]\n",
      " [ 7.93952558e-01  3.50946872e-01  1.18853643e+00  3.85961318e-01\n",
      "   2.96701863e-02]\n",
      " [ 7.26183347e-01  4.41639937e-01  2.71711699e-03  7.33393633e-01\n",
      "   4.55176129e-02]\n",
      " [ 4.89263105e-01  4.20725095e-01  7.56036398e-01  6.24033457e-02\n",
      "  -5.38302416e-04]\n",
      " [ 6.09810836e-01  7.55780427e-01  1.03636918e+00  9.08549910e-01\n",
      "   1.91844947e-03]\n",
      " [ 8.31578328e-01  8.75528332e-05  2.93543168e-01  1.10037225e+00\n",
      "  -2.65884776e-04]\n",
      " [ 4.26650967e-01  5.53761974e-02  6.52855369e-01  6.43132832e-01\n",
      "   1.47569255e-02]]\n",
      "Right factor X:\n",
      "[[ 1.07015116e+00  4.25961964e-01  1.59511553e+00  6.26808607e-01\n",
      "   8.98124301e-01  3.62801718e-01  9.53757673e-01  1.88661317e-01\n",
      "   9.64559055e-01  1.43675625e-01]\n",
      " [ 8.72908811e-01  7.03553498e-01  6.45229205e-01  1.10121868e+00\n",
      "   9.93621271e-01  3.12383803e-01  7.45085312e-01  1.25155585e-01\n",
      "   8.84272390e-01  7.94988511e-01]\n",
      " [ 1.41086863e-04  1.70049131e-01  2.73427259e-01  2.50933223e-02\n",
      "   8.38007474e-03  2.51575697e-01  5.99473425e-01  1.39362252e-01\n",
      "   5.06840502e-01  4.22844259e-01]\n",
      " [ 2.70906925e-01  5.46340550e-01  1.04256418e-02  4.63290841e-01\n",
      "   1.39889787e-01  7.65220031e-03  2.22742919e-01  3.60875098e-02\n",
      "   3.41601146e-01  2.72448408e-02]\n",
      " [ 5.44108256e+00  4.62667224e+00  6.26354249e+00  7.23656013e-01\n",
      "   1.81220987e+00 -2.57729003e-07  2.90739234e+00  2.81123997e+00\n",
      "  -2.15606388e-06  6.43189790e+00]]\n",
      "Residual A - Y * X:\n",
      "[[ 9.02157264e-04  5.23117764e-04 -5.79950842e-04 -5.74317402e-04\n",
      "  -4.61768644e-04 -5.28680186e-05  1.62394448e-04  2.76277321e-04\n",
      "   4.85227596e-04 -5.60481823e-04]\n",
      " [-2.33027425e-04  3.21455250e-04  2.17040399e-04  1.56606195e-04\n",
      "  -2.41256203e-04 -1.01386736e-04  7.36342995e-05 -1.73587325e-05\n",
      "  -5.22429324e-05 -2.04432888e-04]\n",
      " [-8.35846517e-04  2.46121871e-04  5.93720663e-04  5.38806481e-04\n",
      "  -8.42363429e-05 -1.36215640e-04  2.31633730e-06 -1.52108618e-04\n",
      "  -3.23620331e-04 -5.42078084e-06]\n",
      " [ 2.62860853e-04  1.83780003e-05 -3.20542830e-04 -1.49712163e-04\n",
      "  -1.31334078e-04  8.78805144e-05  1.46798183e-04 -2.03546983e-05\n",
      "   4.79256197e-04 -5.81320754e-04]\n",
      " [-6.22557723e-04  6.31892711e-04  4.34719938e-04  4.01388769e-04\n",
      "  -3.52745774e-04 -2.12014739e-04  8.42548761e-05 -4.17321003e-05\n",
      "  -1.50760383e-04 -3.01455643e-04]\n",
      " [-8.46202248e-04  3.61714835e-04  6.15005890e-04  5.85452470e-04\n",
      "  -2.39872783e-04 -1.59000367e-04  6.24749082e-05 -1.69461803e-04\n",
      "  -3.16622183e-04 -8.20910778e-05]\n",
      " [ 1.15561552e-03 -1.28864368e-03 -1.77288000e-03 -5.10264071e-04\n",
      "   6.38713553e-04  7.17730381e-04  2.05892579e-04 -2.69449092e-04\n",
      "   1.71225020e-03 -1.13410340e-03]\n",
      " [ 1.57913703e-04  6.21168134e-04 -4.04695033e-05 -1.48187018e-04\n",
      "  -4.38037868e-04 -1.45409129e-04  1.34145488e-04  1.47289692e-04\n",
      "   1.98184939e-04 -5.09549810e-04]\n",
      " [ 5.51365483e-04 -1.32683206e-03 -1.26345269e-03  6.01647636e-05\n",
      "   9.72529426e-04  6.10472383e-04 -1.48674297e-05 -3.54468161e-04\n",
      "   9.92202367e-04 -1.42249517e-04]\n",
      " [-1.63514531e-03 -1.59800828e-04  1.08957766e-03  1.01954949e-03\n",
      "   3.41048252e-04 -1.06257705e-04 -1.57094132e-04 -3.64204427e-04\n",
      "  -7.26930797e-04  4.63755883e-04]]\n",
      "Residual after 30 iterations: 0.005008539285258121\n"
     ]
    }
   ],
   "source": [
    "#\n",
    "# Plot residuals.\n",
    "#\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Show plot inline in ipython.\n",
    "%matplotlib inline\n",
    "\n",
    "# Set plot properties.\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')\n",
    "font = {'weight' : 'normal',\n",
    "        'size'   : 16}\n",
    "plt.rc('font', **font)\n",
    "\n",
    "# Create the plot.\n",
    "plt.plot(residual)\n",
    "plt.xlabel('Iteration Number')\n",
    "plt.ylabel('Residual Norm')\n",
    "plt.show()\n",
    "\n",
    "#\n",
    "# Print results.\n",
    "#\n",
    "print('Original matrix:')\n",
    "print(A)\n",
    "print('Left factor Y:')\n",
    "print(Y)\n",
    "print('Right factor X:')\n",
    "print(X)\n",
    "print('Residual A - Y * X:')\n",
    "print(A - Y.dot(X))\n",
    "print('Residual after {} iterations: {}'.format(iter_num, prob.value))\n"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}