sectio4_factorizationQR.ipynb

{"cells":[{"cell_type":"markdown","metadata":{"id":"nHyuK0D5ECYp"},"source":["# Factorización QR\n","\n","\n","Algoritmo de factorización QR\n","\n","Para cada  $j$, calculamos una simple proyección  $$v_j = P_ja_j$$  donde $P_j$ es la proyección en el espacio ortogonal del span de $q_1,\\ldots,q_{j-1}$.\n","\n","\n","![Captura1.PNG](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAbEAAADrCAYAAADnsqiUAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAADt7SURBVHhe7d2PU1NXwjfw9z9hB5RRxxaVMnXYKqNjtTKWVaZSf0xFH6poV7DrgnYK2AV8K1AF3j5CbWg31DZakW2wK7hi1/iUtEXHdAt9YAQrItigUFBIESRwv++5NzcxhIRfEkjw+3EyYy75cX9/7zn35Jz/AyIiIj/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/1VCGWe6NfeUw5axcaDaeQl5KIN6OXIzhwOV6LO4D8KjMkqwl5i/6AoIC1yDNZ1DdMtX60nH8HK+asxl/ON2NInTo5Fphy14r5/QNCc02wqlNdSX2/wnROg8ykeLzxSiiCFq1FbNIxlNU/hKS+xmdIvWitLkVe8pt4bdkCdV6PQGf6bcx5tZpyESrWRdCiXJg8rYxJGGr5Gn+JWIwVb3+NFrcbzIo2fYKyHYLi9GhTp067adl/p5Pr/j2Vx44n0/Ed5C98L8T6GqBPfhXB8skmIAihr2zBzv+KRkTgSmQYu6bpJNABQ8py5ftfTK9Cjzp1csYKMQlW8zfIjFyovMaxzFvXipN9EMKzqtGrvtI3DMBcfgBhyrwuQET0dtuFRsCrYnuMvaa8E2JD6DYcwnPy575wGMYed6c1NyHW14SqU/+NtLhMVLRNYaIqJPHxVTh9PA07UyqehOasD7GpPHYEt9toir+D/JqPhVg3ao5vsgXYkj0oqr6LPvXSXurrxsN+8WRaTgIiWB40wlj5PRofDKjTJmuMEBuoR3FMiPL34MhMnG/oUl8j5sFyF413fexEN3QDxVGi9BUQjn36ZnVeh9DXYUaHfWONwlslMaX0bvw3jI329efKTYi16bFTfh6QgLIpDzEPJb9ZH2JTeewIbrfRFH8H+TWfCjGpvRz75sg77CqkGe6LXdUN15PAQAuMxdn460a5NCCmL4rEzox/oN4yqL5hEJaGCuQnbRKlOREUyzYhWfsfcegJUjcayz9EsvJeUarYeBDFNQ/EH9wFjzhwuupQUXAIf7ZXce4uRm3/I7QaP0f2ftvnBwWEIjLufZQ6qgFHC7Eh9BgP40V5vgO3o7hhlDKXfbkX5cBQ8w9kbI1AWO51PGo14rOsvyBGrtZTSnG7kFlaB4vy5fYT6WLsPPMDrhan441l4er83Ue76UtkvrURKyLEOin4Fmbr2CGEwRoU/DFIWc7Yohr1e1x5WldiuiPEMlEm1ltG3KuIeGU3cs7dUD+rFWVxi8Xnv4USUxU++1ssIpao27TLDNOp98VnrsKKjckorLrrWJ8jw9F5u4diXcJRaP62UdkOSqg4To5OD7fVjP1orzmH/ORYRC5aKEqebyGj2IjWPntpz7591yL38g8oTZfnNwtlJWqAOR5iG+hbnfbfNcjSX8Rn6bsQuexVxGd/jUbHPitY76NG/yEObBMlcnn9vZWJz4x3HBd1w5ZvznLEJH0Mgz5TWQfD9tn2n1CWdwCxr4SKfX8j/pz+Ob5r7XXZN+V5v4aKvH1YtyhMrKsPUdHQrb6mb4L7t+fnwx+243dgtP3X4zZy/Q7ZRLaTp2Ulf+RDISahvzobi+Uddf4hGLo91HS7hphcMvjTarzx7hHk5b6P5Ogw8bcgvJh6GV1iz5S6LiPtBXHSDYxE4pFjyE3fh7d0DSI+BtFlSFcCJPiVBGTlfoCMvQeha+wTXzLyIJG6vkOOUuUnDrRIcYAcO4KUtHNolvrQqN2OiK0HkZ17DNlJG20n0xfSYeiST0ruDji7XtQeX6/8LXiXHubRjiT7cs+PxGtRtpKb/HmPG7XYsGwbUrKOIS9rP15bJAeM/SLAHmJBWLwoBPOWRSFaOVmI71uzAdFLwhAZvcY2vwFLsa/813EczE6lZTn4t76HooofYXacLEZbV05hIz8WrcFr8v0/+f+BW6Cp+1282x5iCxEq5jliQ5R68lyItdGvIsz5PXOSUNFuW6OuISaZz+PtJfK6WIAVcWnIsm8X+SGfCB+YoEt/ExHKtJXYmX4UeToTHiqfZif2kepjWCd/f+Aq7PxbOhJfsS3Xiwcuol1ZWfbtG4posS5tF1I5qKw+LQJ6pe37lr2JjNwPoTN1Ou2/8kMEgmP9L8R6zc94LH+k1IHq7Gjls4Ij3kT6+wmIVNbBShysNCvbSGq/iIPyfq0u35G9keo2cbfPurzmhVRUKuvNPu9zlXUd+sp6EQDyZ4rvfbUQNY/kb5ro/u36vBf1uoPYGfemeMTZ7vfK7w+Mg66pH0Oj7b8et5Hrd0xkO422rOSPfCjEnKpfRqtqcg0xsQP/bulznHwlsx7x8s68MBPGXqcr/8h8XO1yrnqw79QLsC73B3QNK4W4HiT9aNLF2U4qG0+gxvmKWZB+t0B8lfqkBWW75ANVvYc34rOcmVGxVw5d578N4qHpjAjkXPVxBqYH4vscyx2EsF2foLrpV9x7IJZbegRLr31+rDDrE5X5DEk3itOH0zpdeRRXewZhvVOKeKW0G4b4kluwSp24mh2lvOa5FIOIqLFJlpuoyNryJBTksNp4FAazvC+Mvq4c22NOAs7eERcM1kbotsmhtQDrtTfExYU9xOZiVfYP6BEXCS2iVDPP6T1Szw/IWTlXvGY5Ug0dwz9X2XceoV4To8zbvF2laJG3rdNyOkpcY1UnDt2E7nX5giEMu/V3bOEh9q/dyvrbimKXCx5bFfgvMN97IEpMY1UnqutfhESTbpctXKK0aBTXAkNNOmyW9+E5+1BmlvfZAbFd99nWgfIaOVi2Kp/tWD7rLZzd5bwvPdkO8+LVCyTpDsri5dfY1/WTebd/jrXpFGKVwIyBpv6RMssT27897++S5UcUbpTXpyjFa+vEkssTR9t/BbfbyOU7JridRltW8j8+FGJON+cnUhITJ31L0/9Al5uKRKUxhPw3+aHu9JYfkB0hn/DEtMAIvPG3z3G1TZ5ncQhXf4AV6uuDl8UivbgabUqYuR6I9rAJQWxJsyMw7STLLVzR5SF171bH1Z2j+miUgxq4j8qkcOVvi7OqxWlH5nTyUx7qcjiW237yUEk9uH1FbsmZ8OQqV34oJ043J9IR68/1Nd1oNHyNMr1+5MPQKJbGbgAPG6tQcmw3VignAnHVq9xkH31duZaYRq4fe4jZ15+794z1GvvfnefBzboYK8Tay5GoLFs8zraqF0BSM85uk0+YI7ev48SrGCvE7Ovfdd4Hxdfut4XathK0qitQai1BrGNe77hZx+PZZwfQWhKvvMY2T272TTfzN7H9281nyqR2GDMilenzYk/h1oB9wUbbf4XxhNgEt9Noy0r+x7fuiYkr0E3KDjv+e2IDTSXYLVcbydWFx7Q4+49c7JzvvNMPoU+ud5fvrSg7ugis13VokjNSbipu/BzpWyNsJw1xwG/W3Rx2hWrb4e0nxYXYpPtl+HwN3MTZ+JeUv0Xu/QDFJaeRF7dUPB/l4HF4ckUdtDIfJucqDcdyuoaY8wH3CE0l+5SWgsGvJCJXewZf5b5puxBQTgKTCTH7ssqvcXnYP2OYfvyq32/7zmEB4mZdCX4VYo6/uzs5hiGx3CwmeNq+kw0xp/e5C7HA/ahovz1i+T2vR/chFry3HO3u5t11/ia8f7tbH1a0V6ba7v0uSUJZi1IGE8baf4XxhNhkt5PbY4r8jU+FmDhSUZEoN52Vd/bhrRNh7UefUm3ivOO1O6qN7CUZqduAtIVOO739feIwtt6vQLJTwFn7+tWD7DHulx90Onhcd/gn1VOuVWSD9Rqskd/3fDauyq0nxRWnITVCvHY8ISbmyVG9txDr0v/5pLWVYzlHCbHBOmhelkuZq5FdLVcEDorS7HsIcSzHZEJsLH24W/szWizqkki/49apt5yqujyvK9n0hJhTdWLs52h4JK5YrKL0kqzeoxpxgtwMTb2bRjWOlphuqqnUezqet6/Ten21CPX2koeb7ei6fPJ9ovXyczfVibYLsCf3Up9LLMd9ZcbsVYX2+XCqchxRnejpYk2ZmWHzN/H92/W52Mdb9Nin3J9ci0OGNmU9KsbcfwW328jlOya7ndwdU+R3fCvE5B3e3W+mlN+JqTvasB2vEy0lu22lqEUbkZx1WG3YIf/ddvIfrC/CBqVxQS5yU7barvqiilD3uFec6LZj3d505ObmIHWrXK1nv7k+cod/0lDA3ljhA6Qf+gcam0uwXSnhheG1pMNPbnyPK8Rk/TBfPmK7Ka28LxSRW+PU3149WQ63B5zUhLOx8slczFP0fmQ7boyL13krxBwnDHk+dzhVATk1OvCwrm7KWTItISZm804JdigXB2J7L9uAra+EIXTZS7bX2JfTcRKV37cWW9MrcU/5NDtRgjBk2qqchzUYEBcc2d8pDYdG276OAFD241hkVP4qz+iI7Thi+aQ2GNJsnzmsYUdgNHKqO8Q6ljDQUIxNyj4jl5AykeWuYUf7v3FIqUof3rAjOPIYqj01OnKZP6llovu3y/OhZpTtlkty4j1yK8v/kht4yI+3kP1v4xj7r+B2G7l+5yS3k7tjivyOj4WYjWS5g6v6QqTt3eHoEeKNve+heFjrLnXHs95FVcHbWLdorq15bsnX0MbJJ1a5auExHtWX4lBcpHqSiMTOdPs9sR7Ul2Rip3ISlk8yu5Dh8Z6YTL73dgXF6W/Z5kk+IBNOovZRH8xVJ7AvUnyO0rz/NC5q5WC1V+O4+yxX4rNbrqJMI0J41xZEKssiwntvKvK036DRopz93Rxwcuh/i8KEKLF8tqbPZy8WYad80lGqorwQYupPGmzzKTdl3oHE1EKcq7nvtGye1pWY32kKMbn00lVTqvwUIVi+F5r+JUw1p21X9Y7lFPNZ/w9kyvuHPI9HDLA1E3HWj7bq02oTc/lnGEnIK7mu7ieyUbav9BD1pe8r+1jwsi3IMoiIdLMdR867PNGMq7oc209HxLzFJOXh7HWz0+fL86VWhSvL9yk079pKmk/mQ6zvtmro1Cbsys9LckthUvZ/2XhO7OIia0L7t8tzx+e5PuRt1zLG/ivPkLtt5G6dT2I7uT2myN/4ZIgR0Vis6Ot7EmlPWifa7wERPRsYYkT+SK7WfX09dqY4/z7SuaqQ6NnAECPyR5Y66LOSlJ447PfckrNLYWrn8UjPFoYYERH5LYYYERH5LYYYERH5LYYYERH5LYYYERH5LYYYERH5LYYYERH5LYYYERH5LYYYERH5LYYYERH5LYbYU5I6avD18UREyD1zz9mEjC+qcNsxCBqNynofNfpcxEd++KT3dn8iWdBS/RXydkUhXu1Nf3LkUQyuoaIoBa8ty1PXRT+adHEIfuEwjD0eRjn3BnmZrp5HUUoMIhw9xE/E08y3BGtHHS5pU53WwwyztqP+khap0THs6d5HMcSmgn1Ih/EMZTLVrE0oe/t9VE55p69WdJi+wmeVt2Afh3dKDTTh0kkNsrbKQ3w4DUHiT6wdaLyUi9echoSZnH50NP4AXdJyp3Uhoa/xAgqLjTA7hhOZBvIyXdMheaHzMCcT8TTzLcK89brLephZkqUVP+mSEMLhWnwWQ2wqPG2ISQ9QU5CLc/LglxOhvG8rgr1ywHfCmL5aHUlYnTTl1DGe/DXEZMrIw08bYjIfWhfq/jy5EHtavrdP2MZ7m4IQm+xx/jRm4junGUNsKjxViNlHdlZHcB63QXRVf4iYJXMRFKPDbb+swWSIPcEQs5mtITbZ4/xpzMR3Tj+G2FQYV4jJI+T+HZkZWchOikHExqMwmPshtX2L/Dh5RN6V2Jl+DCcMLRi74CPB2vJPHM4phmZP2BRXY4rPbv8J5wpS8May7Shu9EplosrXTljystei8tR/Iy0uGbrGh+qoxquwSfMTHqmvGkYJsYWITslCmjKmVxBC1W2rkEce1xxBZvZhJEevRkz2ZQ/VbM7rQp6PGlRo38fOiCjl5ClZmmAs+W8ciF6PFI0GB5Th98MQk/Mt2pWPE+8xV0GTkYkceZj/ZVuRY7hrCyFlhOlMJKZ8hJLPMxH7urrvWRpxSZuFxMgopHz8IeKXBNlK3o+dQkzqQdO5NGXo/+f2fo7/7RDLJabdOvU2lid8iQaLczX2ZOZbvMtSh7MHdmPfkaPITo7Ha2tCnPYJD8eN43Mj8ecj7yM+YoHyua+llqJemScP60PM+21jKQqSN2L5uyegS34Vwco2y4WxXd0J5fVVcgg7EzKRl3UQb0WvwjxHiE3uc8c8zq33YSrag7CABViXLdaN2Ackyw2c3rsKMXnfKc89bl/76NdJqcjOfQ87X9mMtNI69Jjdfaf82lIc2rUP6cc+QHrcBsSo68z9/vA56m7ocSj5A+i+OomcvZvx13NPe8E2tRhiU2HMEJMw0FCMLbv1MMsHrvQLdDEL1ao6qzgPJogdbQJXS311+CxNi5qeWzi7LcQLV8wDaC2JR9DLGtS7vdXWi9rj68U8i2Ue4xG8txzt6rtG8sWSmBVmfSKCl2ZCf+HvKKyoRsXfVmJeUiW61FcMo4TYXKxI+RpNlgFYmr5GasRcddv2olH7Fvbo74g9QGz2Jh02BYRgs+6mmwsV13Uhob86G4udTp6256GILaqBRRpEZ+W7eC5gB3RN4hgcqEfx5mSUmQfEa/txW7cDQYFxyt8G6zVY4/icVpTFLUZIulFsRftnhiCm4BrMDRdx8lwdLK4lMekOyuLDnKqWLWL7JyC7+oH8xMUE51tqRUXiq4gvuaV+1z1UJr/kCHPPx439c8Ows1jMs9SPdrlmInAh1mt+xuNR1od4I7Kf/wPmbfs76uSTd2clDsxfiE26X8TcDsBcfgDhu0rRolxs2OdXXZ5Jf+44jvPHP6Fg1QK8mF6FHvm59Csq/pKJyk5xEHr83j4MNJVgd9QxXO2RD1a1Yc2cd1DZ1T/iO6Wuy0hbugPFDb3Kcww0QhcbhrDE82Idu9sfLuHsX1chtqRZLIN4f4seOQyxWWjMEHuEek0M5kUn4VhuLvJys8RV+WKxE+5HRfvIHW1UUgeq83NxrkWUkJTvnYs1mjpxqLkjrrrutqBDORgn4gGuZkUiPKtanOi8yRdDTF52MU9/jMaBU/Xok0dQjgpVT0RujKhOfIyWkt3iKlycpG+aoHlZlA6SssQ2F9v9aBJem+Mp2EeuC9dqrNGeK0E1JxrJR+X9KxfHkqJF6SEMieVm8UIzTBeq0CiHbMslZEeKUou6r7qtKnMNMbHkj0z5WBWwXryu23ay3foRah+7368mPN/2QFM4r4fRjpvBkfMutaBsV6hy8fXzz6OtD5flc34+WCe22eJh23vE/E7mc8cTYnKA6vdh3pwksXyi9GY+h7/m2I5Bz9/bIC4oorFYHKuOM7G1Azdvdrj5TvUCbdjxpk4LiMfZ1oFhy2rTLbaHuGBdsh05+p/QPtCLhw8fq3/zDQyxqaDurJ5DTL76DfdwIpxISawf5opcHDxyTNmR89LfRMRo92OGbkL3ehQyjJ3qhHFSrghXItXQoU7wFh8MMWXZRUkqphiNA5I4L5Zge+BWz9Wqbu6JOU4Eeq34m/MJejRPE2IPbPuQx3uj4mKmoQL5qe+jqEKP41tCJhhigrIvhWBe/Fe4+cMxvKG94bHae8LzPWzfd14Pox037ubd/t4slJWMsj5GC5vRtudY63m0zx3nca6UlF4IxfaSetzSpaOwRl429b1uv1deR4uHh5iD63e6P96c1+PIdSrmyXIDZanrRdCJC7CI/dDVPXR/QTdDGGJTQd1Z3YbYox9xIl+HL5PCxQlArRZRSV038GNzz7h2bvlquK9Bh3f+3zX02D9DOeBWipByW9E1abarPtvVoHuztzpxqFGL9fYShzgJtJcnYd6qAo+lDts2cD7pqVUycnXOzQtInhOG3Wp1os0gun6uQbMIyOGeJsS60VX5jriC36dWN6mkTvxsakZ3TSGiFr5rq5ZST3oTDjF5vg3peDEwEm9sPYCzLZ6vxicy392GQ3gu4FXxf6UCTXBeD/dR6fG4eeRm3rtgTJerfi/gF4/r4w4GRgubbgNS5wfhpdzrjlBwnl/P63mMzxX/xnecqyWfNZsRu+e0Wn0r1r3H772OiynLETRsH5Xvn/0vau89cvlOe6nL+cLKaX/tclO6Fe/5rfmurbq2Ro+cbeEIWpqDq70ejocZwBCbCurOOiLErPdQnfemOIndUurZwwJeQrymCi2WR+iov4jC1BOo7nps20GVHesRLM138Jub/UOy/IjC3R+q9d42tqoY5x3uCcnSjCu5mxEcW4IW9fOkrmsoOpCKwur7TidVV31o1G6dphaP7kNstPmULD/jdOpB5BvvDQsGS90ZpCUWPLk5r/AwXfoNpqI0HCz83tG4wEa91/BqEeqVkOmAISUCLx35El99+s2wE6nDiBCTT0KbEKOpRZ9yv0ecYJbsgcbYDMtAO+orTyAt/3t0jfispwkxCyTzeby9ZC7CdmlgbOnBQEcdKgsykV/djBr5gmO+LcQkyzXkq9WJLQ8f4uF115OW4HoytpPvy8TIpbHhoeJqQvOtlDzk+f5cuY/kuCf2fDau9tvuT7k/buwn3NXIrpYvOMRbe35AzqotKKjpHmV9dEAaLWyk+zCkrhLb7C18ppQ47PfEbN8z6c+1B9EYx7lYCrGaixAVMPzix/P33ke38TBeDFiIdamnUXXdBFPVlzh86Cs0DYz8zo5OeX0vtu2fyif3ol4Ti6g828Wx67ZS9sv8TEd4Sq0liFUDb6CpHIeTjz9pxDRDGGJPSe6xo0KTpLTeClq0Fm/EvYmdymM7Xlsmt5hSq6LkVmrH92BFoFxCWYAVuwpRpW58qf0KDkcuRFDgBmRetrc4srOio+YfyJavgCKScbKmXezY8rQL0KbIRXyxQ6afRGV9l9NJXSZfla4dVhUz1KLHniVB7kuMdlITzsaGjXKfbYrIP6o1liIzOgRBc2KR969raJRbvgmjzaftbyHYoL3hNH/9aNH/BWGB21Dc4Fzt52H6UDPKdosTpetVsbLsS7HeXlWmPA8d3trQ1cAtVGTEYsXKN5F+7Biy//YuMj75Hm3KfUj5ivhbFOxapVTFBAWuQvzxb920TuwXJ+cLyJN/+D1nBwou30Bb2w0Yju8QJ6BQvJF7AXW3G8T+Iz9fKk7oRrS03VKfhyK2oBrt1j6YqwrVVnpytc8eFFTJ+5Jcgj+DffJ0+ftzS/DVsU1ifsKw8UAhvjgaq35HOX5sFaEi91BxMQ9vzJEbKIgTVIN8b8VOvke1CwcqnS8gnInl7ZjofMsXGqds86fM9w78ecsqrNubDd2VWyL4PR83thPuAkRsfBsZx3KQGhePzPJf1JOz3KrRzfqQeySpKhAn4j/guV1F+OGO2fF83rYTuN4ut3z8GbqE1eo2W409b21CeGQCsr/4H9y29E7+c0c9zp3I+92eIzAM68DAw/LIf1JaU76LSPs6SiiGSXyf8qcR3ynWd30p0qLFuSolB7npyXg7z6Dsk5LlFgy5LvuDHGJ50VjxX5nQlpTii6w0fKC0ihyC5fqHWDdnI/KuT21N0EQxxGar/uvICxcnb9d7Ob1GfHDUc4MNqb0c++ZEi6tZpyvzmTDGfD49EfI5J3yqWsTnySWxHTnT2w3WKEaWGuhZxBCbpZSqxigtGp3PN9I9GD/6+yjFf9s9oJAxqou8bsz5fFpiOY1aFDp+Z0OeiXVVfQr5mi9wJjcBO7T1cLorM6MYYiRjiM1K8n2tWKw+/hPG1Ri2pwoZ4ftQ8h8D8hNHqTqjZ5BaJatUU51Bo490bi1X45dlbbNVfWXpcEGpZqdnEUNsNpLr1LfvfvKDxrEovQ1cwLnK62hRejsgIvIPDLHZZKgZFUc+wcVvinD4ZC0svDQlolmOITarDMLSWodapVUREdHsxxAjIiK/xRAjIiK/xRAjIiK/5XMhJo8tVJq+G4nZHyMv5QhK632rs8kJsd5HjT4X8ZEfjr9vQKkbjecOI2ZREGxjUx3GuYZu3uMiInLDx0JMHqcoWpy8E1D2a73Sa/boHcg+HanjOr4s/gZN3vjty0ATLp3UIEvuSmjcHdwOosv4AbYrfaBdx9WLx7FTHpju1ULUPGKMERG58q0QU8ZukvsGEyFmvgdTUSb+b/ktL/UQMIQeueNM+4B2XuG+g1uP5EHwDn0MkyOw1PGFAtajoMa7I3sREfkj3woxtdfnsYcr8BcTDbFuNDd3Dqs6tHWtEwNNvdvB8YmInmk+FGKdMH3+jjLybVDASuxMz4fO1ClO7L1oNX6O7Izj0J35FDkZBThf16n0eSd11eLL5DXi9fJQGLWoL30Pry1aik2anzDqKV++V3WuEKlbVyu9oXuvO9MJhtgI6lAKo41nNV3EOqutPI2C1N3Yo7uBfvO3KExYj4iYItSyqpOIZoiPl8Tke0RHsMIxaJ46qNuceOhuyTFlH2huLpZvy0bpv88gIzIEYSmXMdbgAMq4OB5LOD4y6KN0B2Xx63Gw0uwTDTsksx7xgWtxRK/HiYILqC1/DyHq2EJERDPBt0NMakHZrtAnz5WXyNVr9pFX7SH2ZCA6qe93/D5ivCZXEnqrcxDu9RFKnybEBtFTnY8dmVdcBm6cKeo6CwhHVPIZNPQ9msbBM4mI3PPtEHNzj8wWYmKaMmDiyBAbH1sryOdSDLCNCestkw0xeTDFShw+8KXP9BpuX2dBgdttHQu7DiBJRDQDZl1JbFwG66B5+SXsK//VQzXdzFYnSpZanDz0KUw99mo6EWq/3cdvY5YwvUhpORqCVWK9yxWwPjN4JhE90/zknpj9ZOnpntjEQmyoUYv1ATu82LTezkOISb/BVJSGg4Xfj6wqtDahLDEOGXojTCaT7VF9Hvk5FTBLYn1c1+LggSJUt7ssrfQQdacy8HZ+1fDPnOB0qesaig6korD6/rCAl5p02BSwGZp6uan/ELoNh/BceCa++kqHC2ZfGSaRiJ41PhRivWj99kO8obROjMA+bRWaHogTtTzW1ZWTT1onpufhrOm+ElhPWieKElF0Ns7XmOEyGL8bj9FSshvBL2tQ7832CNYONBpLkRkdgqA5scj71zU0dqjraqgZZbtfcgprlXQPVZkbEOxU0rM9lqqlRvsAhXJrzFb1TSr1M4P/pEWj83JNcPpQix57lgSp1bV26jpzjBStPl+0BTkcHZmIZpBvlcSmg/yD4oSXxj/qsdd0wZhzYhINSyT0GgtxrPqB+twLeo344Gi1uKwgIvJtz1yIKfdyFu5D2YxWgVnRbtSicBKlGKndiBOFV2D21v0xURo0fvR3GMx+dnFCRM+kZyTEOmFMj8buM9X4Pv8dHL7MKjAiotngGQmxQViafkDFuUu42sJRj4mIZotn754YERHNGgwxIiLyWwwxIiLyWwwxIiLyWwwxnyPB2qLHX1IvjdkTPxHRs44h5mMky48o3BgywQ6NiYieTQwxXyJ1oDp3K8ICF2KT7hf+FICIaAwMMYcBtBmLkBoXidAXDsPYI3cS2IdGbSxeTDwPs9cTpQ8t547hgzMFeCvQTd+IREQ0AkPMha239vUoqJF7DpTH9bqItC2fuHQWPFVDtdhJ6Gs4hUMf/YieFnnE6bXIM3GIEyKisTDEXCnDwaxEhtHWrEIyn8eRk/WinOY9Utf3yM/8J1qsIjSV8dJioKmXh5pxQ7LgbvNvvF9GRCQwxFwN1qDgj6G26jzJjEtHtU6DU3qB9Q4qstOQlZuLvNyjyIhbKUpwLkO0OBkSJcXNS+3VnUREzzaG2AitKIuTh93/GfcufYQT17vcNLCYupGfG4oz8aHjO9RBPhdmwjjhIVqIiJ49DLERzKjYG4aQPe/gg09/hMVrWTIIS40Ge3J+QI/jOx6hXhMzciRohXh9qwF50eHYXtLkJliJiJ49DLERLDDlRiIsvgRNA16KCqkdNaVZiF2yACuSvkBNh0gsMa224lOkRi5E0JxNyPj8G9Q/cKnG7DUiY+EO6Jpm2zonIpochpgrazPK0o7B0O5rTSck9Jvy8FKUFo28HUZEpGCIOZEsdTiblomTdQ99sLpOrmrcgvXaG2CGERHZMMRklmvIiwxD6LJtyDH46KjPQzdQHLUJBTX8/RgRkR1DTCZ1o/HKZZ8e9VlqKcGO14vR6K37dEREfogh5uOGWiqQpSnHvwuP+mg1JxHRzGGI+TrJgtbaOrRavPiDayIiP8UQIyIiv8UQIyIiv8UQs5Meor40E9v35uBk7iG8V1rnxd46vMR6HzX6XMRHfuimxw8P5EYt5w4jZlEQggKCELrxMM41dD/D99760V5zDnm7tiOfIwkQ+TyGmELC45oCrA6Qx/G6iSZdHIID96OifQrvQ0ntMH15CpVNv6sTpthAEy6d1CBra6iHbqvcGUSX8QNsTz2NquvXcfXicexcEoTgVwtR88gXY8yKDtNX+KzyFvrUKVPrEW5fOoW/H9mGeRwOh8gvMMQU8uCXW0VJRA6xZnRd1yLtcDluT2Vz9p4qZLywGJt1N734Y2W5y6y14w8x6VdUHPoYJkdgDcCs3ydO4Pbx1HxNJ4zpqxH8ug5NXvzFt204HIYYkT9giCnUk78SYv48ovJEQ6wbzc2dw6oOxxzP7BnAECPyHwwxDOKh6XMkRy8WITYXEXHvIU9nwkNRXuprNeKzrCMo0J2GNvt95J+vQ5dVnPKlTtSdPogV8nArcaVorC9FWvRShMYUoXZENVw/2q5/hYKUbYj4k7f7PZxgiI0wiK7KdzBvVQFqH09hKfSpSbC2/4RzBSl4Y9l2FDd6pzLRjiFG5D8YYoqRJTGp61tkRizAS7nXRQyJ52Y94gNDsV3XaBvluU2PnXKI/XEHckov4Kv0DQhe8h4MXW7uo0nNOLstFGs0dSIm3Jm68cmeKsSkOyiLX4+DlWYfbNgxgNaSeAS9rEG9l38yxxAj8h8MMYVriFlh1ici2Ll60WpC3iIRJOF5MPWLU7w9xOyBIfXD8ruH5Og2IHX+WmRXP1AneMvThNggeqrzsSPzCtp9L8GEB7iaFYnwrGoR+d7FECPyHwwxhWuIjSyZOUIsIAFlbSIhXEPMIwm91TkIn38Ihm5v9z8/2RCTYDVX4vCBL9HY55MJBjz+CQWrViLV0KFOcKO9HImBthKr50cYEsvN6hvcY4gR+Q+GmGIKSmIeWVB7PBrzEspHKeHMbHWiZKnFyUOfwtRjr6cTofbbffwm3//zEYP1GqyZk4SKaRjnjSFG5D8YYoqRJS/bPbGFWH38JzyWn3u6JzZWYAzWQfPyYmzSVqLi9Hfo9GoueAgx6TeYitJwsPD7kUFqbUJZYhwy9EaYTCbbo/o88nMqYJYGlZ8bHDxQhGrX8JAeou5UBt7Orxr2mZLlZ5xOPYh8471h99U8Tu+6hqIDqSisvj9s+nDqTyBidLjt+UVThiFG5D8YYuLU2dd6BXnyj4QDgvBighbGpgdi6iAsTVdQ7GidmImcEhPaXVsnBm5C9vkatHmqhuu6jNQlcxG27UNUmb24rqwdaDSWIjM6BEFzYpH3r2to7FC/b6gZZbtfEsunVoXaSfdQlblBlDiHl/aCApZiX/mvYh30o0X/F4Q5l0jt1M8MlltcOjW0GGrRY8+SEGzQ3hjWiGX06UEIitOjTZ02gtSEs7FhozSMmSr96Gj4DvqMGLFOQvFG7nlcbezwzfHliEjBEHtmdMGYcwJXeydalJHQayzEMW82Suk14oOjnhtsSO3l2DcnmgOCEtEIDLFnghXtRi0KJzFqtdRuxInCKzB76/6YKA0aP/o7DB5LqWLey5MQEq+HeRqqEonIvzDEyDfJ3XSF70PJfwzITzw6SsgR0bOMIUa+SerBbeMFnKu8jhYOCEpEHjDEiIjIbzHEaMKGWg04kZuLPJdHvr4ez263wUQ0ExhiNEH9uH06BwVVzei5U4Lt899BZdcAeqo/xP6SplF+60VENPUYYjRBvbhpaoBFbTX43C5bq8Ghu7X4aRp60yAicsYQo0m6j8qk1dg+4dJXH1r07yGt8r76nIho8hhiNDlKz/xbJzi21yAsNScQE8gunYhoajDEaBLkXjwyETLBsb2kru+Rt3EpggN2QNfE/YaInh5DbEbJ/TP+D3RZCVi39D0YOu/AkLEB81z6I5xukuUmKnJTkaH5AkXZKfjzhjCs197Ak4FkemDKjXIMGDou1iacy8zHmY/3iBBz6cORiGiSGGIzrhtXs1Zj3n4t/vlJITRHtuDFxPOT7GLp6Yd0kYdl+Sx+F/Kvdyn3umw9usdAU682npc7Gq7SYt8LC7Euq+JJJ8OjsqChOAsnatrRIo/OPOHxzoiI3GOIzbShGyiOCsHyyP34rM7lPpH1NzTftbg0nBClt7st6PBKX4bdqDm+Bauzf0CP8vG2Fojz7GOoTcoguqo/xuFzTeLT1KFiRquGlCy42/zbhPt4JKJnE0NshkktJdgeGIIYTS2GN5EYQo/xMMJf16HJeUDooZvQvR6FDGOnOmHq2HqL3yxKXWp/8gP1KI4JQUi60WMP86OTR4z+F3IOHFZ/EP0edi6bO+qwK0NNOmxeehjGHm+Pgk1EswFDbEapJZ3576Kycypugj1NdaIoMVW+g3nPZ+OqUurqh7n8XbwYsByphg7bSxT96Gg2wzKegllfHT47+LHTiNGtKItb/BShSEQ0HENsRnXAkBKB8CzXsbREUNw4g+SFa5HtNI6XZGnGldzNCI4tQcuU1yaqIaYEai/aqk/iyK4ozJuThArnHzHLvcsvfXPs1oXSA9QUHEROdeeT6lBllOu5CM01uakuHISl1YC86PBJ/PaMiJ5VDLGZ1FuN7KURLiUdG6lJh03zD8HQ7Vyt1gVj+lps0v3ilZO81PUdciIXImhRDNJKr8Ok3TqpwJQ6/oPSrO0IC1iD5JM/okO8X+qoQYU2BesC/4B5GzPwRWU9Hrp+bq8RGQvZ/J6Ixo8h5pNs1YyLUwzoVqco+q8jL3zbBH9gPElKg5PFTk3rJVg7anE6yV3JcSpI6Dfl4aUoLRp5O4yIxokh5pPkasYoJGovwtT8pF/4wXoN1kzTSd7W4MSpab1M+gW6mNVuS45P7xHqNVtcfo9GRDQ6hpgvGhSloD+FIyb7CsyOpvR9aNTGYvXxn/BYneI1yv2srQgOiBWlvichprRefN61inOKKCW/TSioYXdURDR+DDF/ITXh7PbdKG7wdru+33HbeB5ler3tcaFGuaclN/nvNhzC4r2foNJ0BwPKa6eOXPLb8XoxGgcmeAOOiJ5pDDFfN9SMiiOf4OI3RTh8snZ8Tdu9Qi4JbsNzG4/hW/PUbfOhlgpkacrx78KjOFn3kK0SiWhCGGI+T256XofaVteeO2YJyYLW2jq0Wmaws0gi8lsMMSIi8lsMMSIi8lsMMSIi8lsMMSIi8lsMMSIi8lsMMSIi8lsMMR8gtVWhMOVNRC5a4xgnbLBBiw0vHECFeap/VjzN3yf14PaVL5C9NwrhKZfRaf4GmZFh2KC9gYk0qpcsN1GRm4oMzRcoyk7BnzeEsYsqImKI+QylX8IQLD9eYzu5W1tRmbpveN+FwmDNcSx3Mz7YiEfgflS0jxIT4/y+KdFfjeznw/HXL0rx6fECZG1cg7fLWyEpnQq34K7rb8RcRneWLLX4LH4X8q93Kb+Vs5pyERrg3K/jBMY4I6JZhSHmM2xD9zsGjJTu4Nz7X3qxG6bp+76hRi3WB4RjXcLnqOtz/vxOGNOjsFl3c1iJavjozt2oOb4Fq7N/QI/yVnUg0fA8mJTBO4XxjnFGRLMOQ8xnqKMyK0P3ixP1pY9wQi15eMd4v29QhMZ+BLsr7Tk/PJb8HqOlZDeCA7dBUzfxzn2VTofnbBalLrXPyIF6FIsSJEeHJiIZQ8xnWNGmT0BQlBYN977BsY+uqSWP4aasOnGc3/fUpF9RkbAUzyVXotP5860duHEqCSFLc3C11/4H19Gd1dGmn8/GVaXU1Q9z+bt4MWC5OhyMt8c4IyJfxxDzGWqJZ+EOHMjRocbrfQlO0/d1G5A6fy2yqx+oE+z6cVu3A8+5Dvw5bHRnNcTmv4vKzl60VZ/EkV1RmDcnSQS0esfMq2OcEZGvY4j5EKXBwpJ9ONvkhcYVbnj/+yT0VucgfL6bMciUElqUS/iMHN1Z6voOOZELEbQoBmml12HSbkVwbAla1MKbV8c4IyKfxxDzGf34VX8E6YY2N/elvGG6v8+FXEJ7PhHaShOaHY1JxhjdWRk4c7HT3707xhkR+T6GmC+QHqK+5H2kTdd4YdP9fW4ov0ubvxk5396FvSn9WKM7ywNnbg90blrvnTHOiMh/MMRmVBdMuTEIXRSBTdmXYbZ6O1Gm+/smZtTRnaUHqCnYiuCAWBQ3Tk91KxH5PobYjBqEpdGIS1fvTFOJaLq/b3zGHt35d9w2nkeZXm97XKhBh2/lLxHNEIYYzTyO7kxEk8QQIyIiv8UQIyIiv8UQm2LSg0YYNIl4Ue41Y1EuTI6md0RENNUYYt7QpsdOhhgRkdcxxLyBIUZENC0YYt7AECMimhYMMYUEa1cNzmakIFt7BsXZmcgvvwGLZEVX3RkciJiLoIA1yNJfQFHSeoQGRiA2qwJN9rGxpG40ln+IfZGr8FpSDj7T5WArQ4yIyOsYYjKpFRWJyxEUo8NtkUtSkw6bAtYjzyT3r96KsrjFIsSWIl77E7r6xGuTXhLP7cOBDKLLeAQrAuZilTJw4yAspg+xjiFGROR1DDHBNvKwCJ1XUlAs9wihTUGkCKU1mjoRUfYQWytCTe7TzzYiclDAYuzUt4rn91GZFC6er0Z2tTqoCKsTiYimBUNMUIYkkUPnZQ3qR3QaMUaIDdag4I9BTn8XGGJERNOCISZrL0dioAgdZfBFe4pJsFrl/48RYkrP6wvEc5bEiIimG0NMJt2HIXWVCKKFWJdxAY3mX9FUXYrPjfdElI1VndiLRu12BAcswLrj/0EfBtBlzMEqhhgRkdcxxFSS5QbOZe1AhCiRBUfsRt6I1omhiC28js7O6yjcFiqez8WKlHK0yC0UrWZUFyVj3aKleC0pF6dLjtpaJ87ZgYIf7opgIyIib2CIecUQ+h7cR9u9B7C3wicioqnHECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FECMiIr/FEPNp8g+ro9UfVausJuQtSkBZ2/h/Ra10qxWnR5v6fDKfQUTkixhiPo0hRkQ0GoaYT2OIERGNhiHm0xhiRESjYYjJpB7cvvIFsvdGITzlMjrN3yAzMgwbtDcwYmSWaTW+EJPaqlCY8iYiF61BhrFTmTbYoMWGFw6gwjzAECOiWYshZtdfjeznw/HXL0rx6fECZG1cg7fLW/E0XR8O1hzHcrkj4LEegftR0e4uLidQEpN+gS4mBMuP19iC19qKytR90NQ/YogR0azFEFPZRncOx7qEz1HnqddeyYK7zb9h2Knf3bQpM5HqRPm1axGSbkSv/FS6g3Pvf4nGAWmMEBuEpe4U9kXIY6ItwIqEU6izzGz5k4hovBhiisdoKdmN4MBt0NSpozO7MdSkw+alh2HsGVKnuJ82dSYSYr2oPb4eQUpYWdF+6SOcuN6llCRHDbGBGzidcx4tVknk8TXkR4Zhs+4mvLE0RERTjSEmk35FRcJSPJdcic4pHDplWqsTRXC16RMQFKVFw71vcOyja+hRl2XUEOtrQX3z7+of5O97FWs0dTN8L5CIaHwYYrJuA1Lnr0V29QN1gqtBWFoNyIsOx/aSJvU+mbtpU20iITaI9vL9CF64AwdydKhxqhIc9z0x6Q7Kdr+F4galQpKIyOcxxET89FbnIHz+IRi6R6lE6zUiQwSErslped1Nm1ITCTE1rJbsw9mmR+oUm/GFmIRHNRq8lf+kBEdE5OsYYuMiod+Uh5eitGh05Jy7aVNtIiHWj1/1R5BuaBtRKhxPiEld1fi08ArMViYYEfkPhti4PEK9ZgvWa284NXhwN22qjTPEpIeoL3kfaSdrYXGTQWOFmGT5GV/+v/NoGpDfPIieH4340SsNVYiIphZDbDyGbqA4ahMKapxaLrqbNuXGCrEu8fcYhC6KwKbsyx5LUaOGWF89dLteGtbQZF5COdpZICMiP8AQGweppQQ7Xi9WfnNl527a9BuEpdGIS1fvuC2BERHNdgyxUQy1VCBLU45/Fx7FybqHyr0md9OIiGhmMMRGI1nQWluHVuceLNxNIyKiGcEQIyIiv8UQIyIiv8UQ83tyzyE30NjB7UBEzx6GmL9T+n3cBU09u4oiomcPQ8zPKU39N+vQxN8mE9EziCHm1/rRpHsLO7zWATERkW9jiDmRLDdRkZuKDM0XKMpOwZ83hHm5Wyl3BtBmLEJqXCRCX7CPU9aHRm0sXkw8D7NzWg38LzR/2o8y84DyVLLcwpUvspEYGYlUw12YLx/BujnbUNzQp/ydiGi2YYipJEstPovfhXzngSQDYpTh/V1JFjOa3TakUAemdOrCydMjeG852tV3uSM16bApYD0KauR7XRKs5otI2/IJ6h0/T5PwuKYAa4d1ESWhvzobi+fsh+6fWhR8nImYFw6gQg05O3fz73mZiIh8F0NM0Y2a41uwOvsHdRgSK9rLkzAvPA+mfteKOrkK702Ep1ehR53iFUr/hiuRYexSnkrm8zhysl6U0+we4GrWVlHi6lCfy+QS21YE/XE99hX/LJ65427+p2mZiIimGENMkNrLsW/O5ict/AbqURwTgpB0oyhbzZDBGhT8MdTW+a9kxqWjWph6nHoJ6alCxpoPcLXXKWSlJpyNXYzgjZ+grs81fImIZh+GGAbRVfkO5j2fjatKqasf5vJ38WLAcpdSjiB3OXUlD68F7sbZlsfqRGdTV50ItKIsbinWa3/GvUsf4YRazWkzhG7De4jIqh4WsrYwfgkHKu85vdaJu/kfc5mIiHwXQ8weYvPfRWVnL9qqT+LIrijMm5OEinanMbsUEnqNmQiJ0eG21ws6ZlTsDUPInnfwwac/Du+lXvlt2HaXYWDkYDuE55bmDC+dDeNu/qdzmYiIphZDTJC6vkNO5EIELYpBWul1mLRbERxbgpYRJ/UemHI3TFOLRXkssUiExZeog1U+IZn12P2nItQPmy7fI1uL51IM6FanjORu/qdzmYiIphZDzJUy2OVi9yf1wTpoXo5FceM0NFm3NqMs7RgMI0qDj9FSsg9bdDcnHjru5n86l4mIaIoxxFzIPWBsD3TftH6oUYv1qwpQ+9i79W6SpQ5n0zLdj1cmh+yfEiZ1/8rd/E/XMhEReQNDzJn0ADUFWxEcIJdMXEPMVgLarHVu5j7FLNeQFxmG0GXbkGO4C9cymJhBN78NGy938z8Ny0RE5EUMMYffcdt4HmV6ve1xoQYdSlD0o6X8ODQXy1GYeRp13hwMU+pG45XLuNpiGVkCU0iwdvyCn1s9/d0dd/M/jctERORFDLFxkCytqK1tHd5C0I+4m39/XyYiIhlDjIiI/BZDjIiI/NbsDDGpC40GDRJfCEJQwFrkmZx/FExERLPFLC6Jyd02LWaIERHNYgwxIiLyWwwxIiLyWz4WYgPoqilFZnI2ir/SIiftv1HR0A1J6kTd6YNYIfcAvygTZVV/R3JkKIKX7UB2+S/quFmDsDRUID9hPSKi9yOnuBi5W0IYYkREs5hPhZg88OPbS0KwSfcLJPTjtm4Hglbmw/RIAtr02CmH2Pw9KP7f39B3vwLJ8+Xnh2DoHoLU9S0yI+aK1x/F1Z5BSJZryI9cwBAjIprFfCjE1FGJAxYgMuVTlOn/geKUKPFc7cfQHmKLcmGS+2NSRj6Wx+ZKQFlbv204FfH3xVnVIv5krE4kIprtfCjE5KFH1orQmYs1mjqM6Ahp1BDrdgxGGZprUvscZIgREc12PhRig2gv349gEUTPJVei094dkjQI66BTdaLbELOopTiWxIiIniW+dU+s6zLS5B8oB25A5sV6mM2/oPr0GXwnj6k1aohZMdBQjE2B4nlkIWr7hmDtqkL2yrkMMSKiWcynQszWwvBrZG+LECWyBVixK29k68Q5O1D40310/qRB7Bw5xNYitbwJfaL81VYtt1oMQ2h0EvJO6ZCntE4MRexxI8x97OmWiGi28bEQm1pS3wPca7uPh6JkRkREs8+sDjEiIprdGGJEROS3nirEiIiIZhJDjIiI/BZDjIiI/BZDjIiI/BZDjIiI/BZDjIiI/BZDjIiI/BZDjIiI/BZDjIiI/BZDjIiI/BZDjIiI/BZDjIiI/BZDjIiI/BZDjIiI/BZDjIiI/BTw/wGkC5TztjX3BgAAAABJRU5ErkJggg==)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"yg3vX5nSF-_I"},"outputs":[],"source":["import numpy as np\n","\n","np.set_printoptions(suppress=True, precision=4)  # precision work"]},{"cell_type":"markdown","metadata":{"id":"TNBGNGYDGS3U"},"source":["## Implementación del método QR básico"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"5EAJlC-TGXzr"},"outputs":[],"source":["def cgs(A):\n","    m, n = A.shape\n","    Q = np.zeros([m,n], dtype=np.float64)\n","    R = np.zeros([n,n], dtype=np.float64)\n","    for j in range(n):\n","        v = A[:,j]\n","        for i in range(j):\n","            R[i,j] = np.dot(Q[:,i], A[:,j])\n","            v = v - (R[i,j] * Q[:,i])\n","        R[j,j] = np.linalg.norm(v)\n","        Q[:, j] = v / R[j,j]\n","    return Q, R"]},{"cell_type":"code","execution_count":null,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"ltdB3-FoGdDt","outputId":"6054c263-b90c-481a-a453-d88e4f9b9ce7"},"outputs":[{"data":{"text/plain":["True"]},"execution_count":3,"metadata":{},"output_type":"execute_result"}],"source":["## creation of example matrix\n","n = 6\n","A = np.random.rand(n,n)\n","\n","Q, R = cgs(A)\n","\n","## comprobation \n","## Con este c[odigo comprobamos si Q*R es igual a A\n","\n","np.allclose(A, Q @ R)"]},{"cell_type":"markdown","metadata":{"id":"7aCj3OZqHFXE"},"source":["\n","El método de Gram-Schmidt, es inestable, pues para cada $j$ calcula una simple proyección $$v_j = P_ja_j$$ donde $P_j$ proyecta ene le espacio ortogonal del span de  $q_1,\\ldots,q_{j-1}$.\n","\n","si modificamos  Gram-Schmidt: para cada  $j$, calculamos $j-1$ proyectando $$P_j = P_{\\perp q_{j-1}\\cdots\\perp q_{2}\\perp q_{1}}$$"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"H8d3lqY3INCG"},"outputs":[],"source":["def mgs(A):\n","    V = A.copy()\n","    m, n = A.shape\n","    Q = np.zeros([m,n], dtype=np.float64)\n","    R = np.zeros([n,n], dtype=np.float64)\n","    for i in range(n):\n","        R[i,i] = np.linalg.norm(V[:,i])\n","        Q[:,i] = V[:,i] / R[i,i]\n","        for j in range(i, n):\n","            R[i,j] = np.dot(Q[:,i],V[:,j])\n","            V[:,j] = V[:,j] - R[i,j]*Q[:,i]\n","    return Q, R"]},{"cell_type":"markdown","metadata":{"id":"Nbt5K0ey1a7S"},"source":["# Nota\n","Si $A_{n\\times n}$, donde cada una de sus columnas son linealmente idependientes  entonces las columnas  de  $Q$ forman una base para $\\mathbb{R}^n $. Además la matriz $Q$ es ortonormal, es decir $QQ^T=Q^TQ=I$.\n","\n","## Pregunta\n","¿Qué podemos decir de las filas de la matriz de $Q$ si $A_{n\\times m}$ con $n>m$? ¿ De los productos $Q^tQ$ y $QQ^T$?"]},{"cell_type":"markdown","metadata":{"id":"JRKrtoSF1a7T"},"source":["# Existe alguna aplicación para QR. \n","\n","La factorización $QR$ tiene muchas aplicaciones, las más relevantes están en torno que puede ayudar a solucionar adecuadamente problemas de mínimos cuadrados, para introducir este concepto, trabajaremos con la interpolación de un polinomio de grado 7.\n","\n","Sea los puntos $x_0=2.0$, $x_1=1.2$,...,$x_{10}=4.0$, puntos igualmente espaciados entre el intervalo $[2,4]$ y sea el conjunto $y_i=1+x_i+x_1^2+\\cdots+x_i^{10}$ para $0\\leq i\\leq 10$. Obviamente el polinomio que que une a estos puntos es $P(x)=1+x+x^2+\\cdots+x^{10}$. supongamos que no concemos el polinomio $P(x)$, y lo queremos encontrar, de esta forma podemos decir que $$P(x)=c_0+c_1x+\\cdots+c_{10}x^{10},$$ de esta forma podemos plantear el sistemas de ecuaciones \n","\n","$$\\begin{pmatrix}\n","1&x_0&x_0^2 &  & x^{10}_{0}\\\\\n","\\vdots&\\vdots & \\vdots & & \\vdots\\\\\n","1&x_{10}&x_{10}^2 &  & x^{10}_{10}\n","\\end{pmatrix}\\begin{pmatrix}c_0\\\\\\vdots\\\\c_{10}\\end{pmatrix}=\\begin{pmatrix}y_0\\\\\\vdots\\\\y_{10}\\end{pmatrix}$$"]},{"cell_type":"code","source":["!pip install numpy\n","!pip install scipy "],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"7k_XiYmYDrXh","executionInfo":{"status":"ok","timestamp":1646148811737,"user_tz":300,"elapsed":6298,"user":{"displayName":"Abel Alvarez Bustos","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"04726289280755598660"}},"outputId":"c7b4e6e5-f6c1-479c-e9b2-8ecdbac6430c"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["Requirement already satisfied: numpy in /usr/local/lib/python3.7/dist-packages (1.21.5)\n","Requirement already satisfied: scipy in /usr/local/lib/python3.7/dist-packages (1.4.1)\n","Requirement already satisfied: numpy>=1.13.3 in /usr/local/lib/python3.7/dist-packages (from scipy) (1.21.5)\n"]}]},{"cell_type":"code","execution_count":null,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"IGxwcbM31a7T","executionInfo":{"status":"ok","timestamp":1646148933286,"user_tz":300,"elapsed":184,"user":{"displayName":"Abel Alvarez Bustos","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"04726289280755598660"}},"outputId":"4fe69f78-1884-4b35-f257-a49deff906ce"},"outputs":[{"output_type":"stream","name":"stdout","text":["[[1.00000000e+00 2.00000000e+00 4.00000000e+00 8.00000000e+00\n","  1.60000000e+01 3.20000000e+01 6.40000000e+01 1.28000000e+02\n","  2.56000000e+02 5.12000000e+02 1.02400000e+03]\n"," [1.00000000e+00 2.20000000e+00 4.84000000e+00 1.06480000e+01\n","  2.34256000e+01 5.15363200e+01 1.13379904e+02 2.49435789e+02\n","  5.48758735e+02 1.20726922e+03 2.65599228e+03]\n"," [1.00000000e+00 2.40000000e+00 5.76000000e+00 1.38240000e+01\n","  3.31776000e+01 7.96262400e+01 1.91102976e+02 4.58647142e+02\n","  1.10075314e+03 2.64180754e+03 6.34033810e+03]\n"," [1.00000000e+00 2.60000000e+00 6.76000000e+00 1.75760000e+01\n","  4.56976000e+01 1.18813760e+02 3.08915776e+02 8.03181018e+02\n","  2.08827065e+03 5.42950368e+03 1.41167096e+04]\n"," [1.00000000e+00 2.80000000e+00 7.84000000e+00 2.19520000e+01\n","  6.14656000e+01 1.72103680e+02 4.81890304e+02 1.34929285e+03\n","  3.77801998e+03 1.05784560e+04 2.96196767e+04]\n"," [1.00000000e+00 3.00000000e+00 9.00000000e+00 2.70000000e+01\n","  8.10000000e+01 2.43000000e+02 7.29000000e+02 2.18700000e+03\n","  6.56100000e+03 1.96830000e+04 5.90490000e+04]\n"," [1.00000000e+00 3.20000000e+00 1.02400000e+01 3.27680000e+01\n","  1.04857600e+02 3.35544320e+02 1.07374182e+03 3.43597384e+03\n","  1.09951163e+04 3.51843721e+04 1.12589991e+05]\n"," [1.00000000e+00 3.40000000e+00 1.15600000e+01 3.93040000e+01\n","  1.33633600e+02 4.54354240e+02 1.54480442e+03 5.25233501e+03\n","  1.78579390e+04 6.07169928e+04 2.06437775e+05]\n"," [1.00000000e+00 3.60000000e+00 1.29600000e+01 4.66560000e+01\n","  1.67961600e+02 6.04661760e+02 2.17678234e+03 7.83641641e+03\n","  2.82110991e+04 1.01559957e+05 3.65615844e+05]\n"," [1.00000000e+00 3.80000000e+00 1.44400000e+01 5.48720000e+01\n","  2.08513600e+02 7.92351680e+02 3.01093638e+03 1.14415583e+04\n","  4.34779214e+04 1.65216101e+05 6.27821185e+05]\n"," [1.00000000e+00 4.00000000e+00 1.60000000e+01 6.40000000e+01\n","  2.56000000e+02 1.02400000e+03 4.09600000e+03 1.63840000e+04\n","  6.55360000e+04 2.62144000e+05 1.04857600e+06]]\n","[   2047.            4868.48584509   10868.43673692   22939.02804367\n","   46074.49704151   88573.          163766.80463165  292453.0984918\n","  506236.93785487  852042.67936983 1398101.        ]\n"]}],"source":["\n","import numpy as np\n","import scipy.linalg as la\n","\n","n=11\n","A = np.zeros([n,n], dtype=np.float64)\n","c = np.zeros([n,1], dtype=np.float64)\n","y = np.zeros([n,1], dtype=np.float64)\n","\n","for j in range(n):\n","    xi=2.0+0.2*j\n","    A[j,:]=np.array([xi**i for i in range(n)])\n","y=np.sum(A,axis=1)    \n","print(A)\n","print(y)\n"]},{"cell_type":"code","source":["P, L, U = la.lu(A)\n","b = la.solve_triangular(L , P@y, lower=True)\n","c = la.solve_triangular(U , b, lower=False)\n","print(c)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"X5-cApI9D99v","executionInfo":{"status":"ok","timestamp":1646148937546,"user_tz":300,"elapsed":250,"user":{"displayName":"Abel Alvarez Bustos","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"04726289280755598660"}},"outputId":"6a22efd5-2eba-41b7-cd7b-062ae92d4962"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["[ 2.76220576e+13 -9.75165846e+13  1.54052241e+14 -1.43411981e+14\n","  8.71287819e+13 -3.60985005e+13  1.03295611e+13 -2.01588641e+12\n","  2.56796792e+11 -1.92824556e+10  6.48134011e+08]\n"]}]},{"cell_type":"code","source":["#Q, R = mgs(A)\n","Q, R = la.qr(A)\n","f = la.solve_triangular(R , Q.T@y, lower=False)\n","print(f)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"DdmcRBDIFP_X","executionInfo":{"status":"ok","timestamp":1646148942122,"user_tz":300,"elapsed":231,"user":{"displayName":"Abel Alvarez Bustos","photoUrl":"https://lh3.googleusercontent.com/a/default-user=s64","userId":"04726289280755598660"}},"outputId":"d012b0d1-665a-48ab-d2a9-08686ed8f89f"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["[1.00105674 0.99623443 1.00600109 0.99436758 1.00344772 0.99856181\n"," 1.00041405 0.99991876 1.0000104  0.99999922 1.00000003]\n"]}]},{"cell_type":"code","source":["\n","\n","\n","b = la.solve_triangular(L , P@y, lower=True)\n","c = la.solve_triangular(U , b, lower=False)\n","print(c)\n","\n","#Q, R = la.qr(A)\n","#Q, R = cgs(A)\n","Q, R = mgs(A)\n","\n","f = la.solve_triangular(R , Q.T@y, lower=False)\n","print(f)\n","\n","#np.allclose(A, P @ L@ U)"],"metadata":{"id":"R2EkWVJODa7w"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"ywCZyNFbIcPa"},"source":["# Householder reflections "]},{"cell_type":"markdown","metadata":{"id":"t2tuhJ4t1a7U"},"source":["## Teorema \n","Asumimos que $w$ y $x$ tienen la misma longitud euclidiana entonces $w-x$ y $w+x$ son perpendiculares.\n","\n","### Matriz de proyección \n","Una matriz se llama de proyección si $A^2=A$\n","$$P=\\frac{vv^t}{v^tv}$$\n","\n","Propiedades de $P$\n","    * $P$ es simétrico \n","    * $Pv=v$\n","    * $P^2=P$\n","Dado un vector $u$, $Pu$ es la proyección de $u$ en $v$\n","\n","Si nosotros substraemos 2 veces $Px$ desde $x$, podremos obtener $w$. Para verificar esto \n","\n","sea $$H=I-2P$$\n","$$Hx=w-v-\\frac{2vv^T}{v^Tv}x=w$$\n","\n","<img src=\"imagenes\\hause.png\" alt=\"term-document matrix\" style=\"width: 30%\"/>\n","\n","La matriz $H$ es llamado el reflector de  Householder\n","\n","    * $H$ es simétrico\n","    * $H^TH=HH=I$, e.d ortogonal\n","    \n","## Reflector de Householder\n","\n","Sea $x$ y $w$ dos vectores de la misma longitud $||x||_2=||w||_2$ y definimos $v=w-x$. entonces $$H=I-2\\frac{vv^T}{v^Tv}$$ es una matriz simetrica ortogonal y $Hx=w$\n","\n","## Ejercicio \n","\n","Sea $x=[3,4]$ y $w=[5,0]$, encontrar el reflector de Hauseholder que satisface $$Hx=w$$\n","\n","\n","### Algortimo de Hoseholder \n","\n","\n","<img src=\"imagenes\\hase1.png\" alt=\"term-document matrix\" style=\"width: 50%\"/>\n","\n","<img src=\"imagenes\\hase2.png\" alt=\"term-document matrix\" style=\"width: 50%\"/>\n","\n","<img src=\"imagenes\\hase3.png\" alt=\"term-document matrix\" style=\"width: 50%\"/>\n"]},{"cell_type":"markdown","metadata":{"id":"bLsZGRZCIxOc"},"source":["## algoritmo"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"uxuEePLIIzNR"},"outputs":[],"source":["def householder_lots(A):\n","    m, n = A.shape\n","    R = np.copy(A)\n","    V = []\n","    Fs = []\n","    for k in range(n):\n","        v = np.copy(R[k:,k])\n","        v = np.reshape(v, (n-k, 1))\n","        v[0] += np.sign(v[0]) * np.linalg.norm(v)\n","        v /= np.linalg.norm(v)\n","        R[k:,k:] = R[k:,k:] - 2*np.matmul(v, np.matmul(v.T, R[k:,k:]))\n","        V.append(v)\n","        F = np.eye(n-k) - 2 * np.matmul(v, v.T)/np.matmul(v.T, v)\n","        Fs.append(F)\n","    return R, V, Fs"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"wbZ0x0pK1a7W"},"outputs":[],"source":[""]},{"cell_type":"markdown","metadata":{"id":"-MQs2Wp61a7W"},"source":["\\begin{array}{ l | l | c }\n","\\hline\n","Gram-Schmidt & Triangular\\, Orthogonalization & A R_1 R_2 \\cdots R_n = Q  \\\\\n","Householder  & Orthogonal\\, Triangularization & Q_n \\cdots Q_2 Q_1 A = R  \\\\\n","\\hline\n","\\end{array}"]},{"cell_type":"markdown","metadata":{"id":"Ct3akZiQ1a7W"},"source":["## Ejemplo\n","\n","Tomado de Trefethen and Bau\n","\n","We want to construct a square matrix A with random singular vectors and widely varying singular values spaced by factors of 2 between $2^{-1}$ and $2^{-(n+1)}$"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"BeZ7cN441a7W"},"outputs":[],"source":["import matplotlib.pyplot as plt\n","from matplotlib import rcParams\n","%matplotlib inline"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"vH4Xc8Rw1a7X"},"outputs":[],"source":["n = 100\n","U, X = np.linalg.qr(np.random.randn(n,n))   # set U to a random orthogonal matrix\n","V, X = np.linalg.qr(np.random.randn(n,n))   # set V to a random orthogonal matrix\n","S = np.diag(np.power(2,np.arange(-1,-(n+1),-1), dtype=float))  # Set S to a diagonal matrix w/ exp \n","                                                              # values between 2^-1 and 2^-(n+1)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"uqF4rCAl1a7X"},"outputs":[],"source":["A = np.matmul(U,np.matmul(S,V))"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"BHwNNJyM1a7X"},"outputs":[],"source":["QC, RC = cgs(A)\n","QM, RM = mgs(A)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"YV06V3xt1a7Y"},"outputs":[],"source":["plt.figure(figsize=(10,10))\n","plt.semilogy(np.diag(S), 'r.', basey=2, label=\"True Singular Values\")\n","plt.semilogy(np.diag(RM), 'go', basey=2, label=\"Modified Gram-Shmidt\")\n","plt.semilogy(np.diag(RC), 'bx', basey=2, label=\"Classic Gram-Shmidt\")\n","plt.legend()\n","rcParams.update({'font.size': 18})"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"8TIFbOK_1a7Y"},"outputs":[],"source":["type(A[0,0]), type(RC[0,0]), type(S[0,0])"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"HyEUCw-A1a7Y"},"outputs":[],"source":["eps = np.finfo(np.float64).eps; eps"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"geZwu30q1a7Y"},"outputs":[],"source":["np.log2(eps), np.log2(np.sqrt(eps))"]},{"cell_type":"markdown","metadata":{"id":"I2TXDrJN1a7Y"},"source":["## perdidad de ortogonalidad\n","\n","Tomado de Trefethen"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"FziIZ_OA1a7Z"},"outputs":[],"source":[""]},{"cell_type":"code","execution_count":null,"metadata":{"id":"5hgATIO_1a7Z"},"outputs":[],"source":["A = np.array([[0.70000, 0.70711], [0.70001, 0.70711]])\n","Q1, R1 = mgs(A) ## Gram smith \n","R2, V, F = householder_lots(A)  #Hauseholder\n","Q2T = np.matmul(block_diag(np.eye(1), F[1]), F[0])"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"g6eYTHD81a7Z"},"outputs":[],"source":["Q3, R3 = np.linalg.qr(A)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"YU8A16VJ1a7Z"},"outputs":[],"source":["np.matmul(Q1, R1)\n","np.matmul(Q2T.T, R2)\n","np.matmul(Q3, R3)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"utaH3mLS1a7a"},"outputs":[],"source":["np.linalg.norm(np.matmul(Q1.T, Q1) - np.eye(2))  # Modified Gram-Schmidt"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"HmM32Kk81a7a"},"outputs":[],"source":["np.linalg.norm(np.matmul(Q2T.T, Q2T) - np.eye(2))  # Our implementation of Householder"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"MQGKHtFv1a7a","outputId":"1db39de0-2f66-47c0-eeba-b6db0342812b"},"outputs":[{"ename":"NameError","evalue":"name 'Q3' is not defined","output_type":"error","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)","\u001b[0;32m/tmp/ipykernel_49082/3621899155.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnorm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatmul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mQ3\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQ3\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meye\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# Numpy (which uses Householder)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m","\u001b[0;31mNameError\u001b[0m: name 'Q3' is not defined"]}],"source":["np.linalg.norm(np.matmul(Q3.T, Q3) - np.eye(2))  # Numpy (which uses Householder)"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"LmAq3U5V1a7a"},"outputs":[],"source":[""]}],"metadata":{"colab":{"name":"sectio4_factorizationQR.ipynb","provenance":[]},"hide_input":false,"kernelspec":{"display_name":"Python 3 (ipykernel)","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.8.11"},"latex_envs":{"LaTeX_envs_menu_present":true,"autoclose":false,"autocomplete":true,"bibliofile":"biblio.bib","cite_by":"apalike","current_citInitial":1,"eqLabelWithNumbers":true,"eqNumInitial":1,"hotkeys":{"equation":"Ctrl-E","itemize":"Ctrl-I"},"labels_anchors":false,"latex_user_defs":false,"report_style_numbering":false,"user_envs_cfg":false}},"nbformat":4,"nbformat_minor":0}