{"id":450,"date":"2016-03-23T15:02:15","date_gmt":"2016-03-23T15:02:15","guid":{"rendered":"http:\/\/drsfenner.org\/blog\/?p=450"},"modified":"2016-03-23T15:03:04","modified_gmt":"2016-03-23T15:03:04","slug":"svd-computation-capstone","status":"publish","type":"post","link":"https:\/\/drsfenner.org\/blog\/2016\/03\/svd-computation-capstone\/","title":{"rendered":"SVD Computation Capstone"},"content":{"rendered":"<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>At long last, I&#8217;ve gotten to my original goal: writing a relatively easy to understand version of your plain, vanilla SVD computation. I&#8217;m going to spend just a few minutes glossing (literally) over the theory behind the code and then dive into some implementation that builds on the previous posts. Thanks for reading!<!--more--><\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3 id=\"words-of-thanks-and-caution\">Words of Thanks and Caution<\/h3>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>I had really hoped to minimize the hand-waving that I&#8217;m about to do. To really get into this, you need to start by reading pages 390 to 450 of Golub &amp; VanLoan (3rd) &#8211; definitely section 8.1-8.3 before starting on 8.6. And, even if you do that, you&#8217;ll probably be at a loss as to what the heck it all means &#8212; if you don&#8217;t bail for beer-thirty, first. Atkinson, Intro. to Numerical Analysis, 2nd, pages 587-633 was also very helpful in understanding the theory behind the code (and it was a trip down memory lane to my numerical analysis course, fifteen+ years ago).<\/p>\n<p>I found a few resources online that smoothed out some of (my?) rough edges:<\/p>\n<ul>\n<li><a href=\"http:\/\/numerical-renaissance.com\/NR.pdf\">Numerical Renaissance<\/a> (Thomas Bewley): pages 108-110<\/li>\n<li><a href=\"http:\/\/cavern.uark.edu\/~arnold\/SciComp\/SC.pdf\">Pages From a Scientific Computation Text<\/a> (Mark Arnold): the chapter on QR Iterations<\/li>\n<li><a href=\"http:\/\/web.stanford.edu\/class\/cme335\/lecture6.pdf\">Jim Lambers Lecture Notes<\/a>.<\/li>\n<\/ul>\n<p>I&#8217;m saying nothing about the applications (or general description) of SVD in this post. There are 100s (or more) of descriptions out there. Application wise, I&#8217;m most interested in SVD for regression problems, but there are many other uses. This post is interested in calculating the SVD. And that requires a bit of theory.<\/p>\n<p>I wanted to actually demonstrate (you know, some running programs) the similarities between a &quot;long winded&quot; (non-specialized QR &#8211; see below) and the &quot;short and sweet&quot; SVD calculations. I may do that in a later post. But for now, we&#8217;ll have to be satisfied with a verbal explanation (aka, hand-waving). A warning: the following paragraphs summarize some 40-50 pages of mathematical textbook presentation. I probably made a couple mistakes and left an important concept by the wayside. However the code that follows should be reasonably clear and is (somewhat) tested. I just wanted to present the mathematical theory with a minimum of, ahem, theory.<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3 id=\"overview-of-qr-based-svd-calculations\">Overview of QR Based SVD Calculations<\/h3>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>Here goes:<\/p>\n<p>It is fairly easy (ha!) to solve the <em>symmetric eigenvalue problem<\/em>. Given a symmetric matrix, <span class=\"math inline\">\\(A=A^T\\)<\/span>, we can diagonalize it with a orthogonal similarity transformation <span class=\"math inline\">\\(Q^T_DAQ_D=D\\)<\/span> with <span class=\"math inline\">\\(D\\)<\/span> diagonal. This is usually done in two steps: (1) <span class=\"math inline\">\\(A\\rightarrow Q_{Tridi}AQ^T_{Tridi}=T\\)<\/span> with <span class=\"math inline\">\\(T\\)<\/span> symmetric and tridiagonal, followed by (2) applying the <span class=\"math inline\">\\(QR\\)<\/span>-iteration to <span class=\"math inline\">\\(T\\)<\/span>. (1) is done with series of Householder reflections. What about (2)?<\/p>\n<p>The <span class=\"math inline\">\\(QR\\)<\/span>-iteration works like this. Using a QR decomposition, let <span class=\"math inline\">\\(T=T_1 = Q_{q_1} R_{r_1}\\)<\/span> and then construct <span class=\"math inline\">\\(T_2=R_{r_1} Q_{q_1}\\)<\/span> by multiplying in the reverse order. We can eliminate <span class=\"math inline\">\\(R_{r_1}\\)<\/span> from that definition of <span class=\"math inline\">\\(T_2\\)<\/span>: <span class=\"math inline\">\\(T_2 = Q^T_{q_1} T_1 Q_{q_1}\\)<\/span> and even better, that is a similarity transformation from <span class=\"math inline\">\\(T_1 \\rightarrow T_2\\)<\/span>. This means that they share eigenvalues. If we repeat this process many times, <span class=\"math inline\">\\(T_i\\)<\/span> will converge to a diagonal matrix since <span class=\"math inline\">\\(A\\)<\/span> was symmetric.<\/p>\n<p>A side note: in general, for a not-necessarily-symmetric <span class=\"math inline\">\\(A\\)<\/span>, the initial reduction is to a Hessenberg matrix (non-zero below the first subdiagonal; an upper triangular matrix extended to the first subdiagonal) and the <span class=\"math inline\">\\(QR\\)<\/span>-iteration converges to an upper triangular matrix.<\/p>\n<p>The summary of this is that if <span class=\"math inline\">\\(Q_{qr}=\\prod_i Q_{r_i}\\)<\/span>, then we have <span class=\"math inline\">\\(Q^T_{Tridi}Q^T_{qr} A Q_{qr} Q_{Tridi} = D\\)<\/span> is diagonal. So, the original <span class=\"math inline\">\\(Q\\)<\/span> we sought is <span class=\"math inline\">\\(Q_{qr}Q_{Tridi}\\)<\/span>. Oh yeah, since this is all one big similarity transformation, the eigenvalues of <span class=\"math inline\">\\(A\\)<\/span> and <span class=\"math inline\">\\(D\\)<\/span> are the same: <span class=\"math inline\">\\(\\mathrm{eig}(A)=\\mathrm{eig}(D)\\)<\/span>. So, <em>diagonalization of a symmetric matrix gets us the eigenvalues of the matrix<\/em>.<\/p>\n<p>One other note: the tridiagonal matrix that is input to the <span class=\"math inline\">\\(QR\\)<\/span>-iteration needs to be non-zero on its sub-diagonal (the diagonal below the main diagonal) elements. If any of them are zero, the <span class=\"math inline\">\\(QR\\)<\/span>-iteration will not necessarily converge. However, those zeros allow us to solve sub-problems above and below the row with the zero off-diagonal entry. In that case, the eigenvalues of the whole matrix are equal to the union of the eigenvalues of the sub-problems. Easy enough.<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3 id=\"wait.-why-do-we-care\">Wait. Why do we Care?<\/h3>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>Well, for a matrix <span class=\"math inline\">\\(A\\)<\/span>: <span class=\"math inline\">\\(\\mathrm{singular\\_values}(A)^2=\\mathrm{eig}(A^TA)=\\mathrm{eig}(AA^T)\\)<\/span> (if <span class=\"math inline\">\\(\\mathrm{eig}\\)<\/span> returns non-zero eigenvalues). If we have the SVD of <span class=\"math inline\">\\(A\\)<\/span>, we can see that <span class=\"math inline\">\\(A = UDV^T\\)<\/span> gives <span class=\"math inline\">\\(A^TA=(UDV^T)^T UDV^T= V D^T U^T U D V^T = V D^T D V^T = V D^2 V^T\\)<\/span>. And we have the squares of <span class=\"math inline\">\\(D\\)<\/span>&#8216;s values in the last equality.<\/p>\n<p>Fortunately for us, <span class=\"math inline\">\\(A^TA\\)<\/span> is symmetric. Unfortunately, we can&#8217;t just apply the above eigenvalue computation to <span class=\"math inline\">\\(A^TA\\)<\/span> because we run into the same condition number (numerical stability) and memory use problems we would hit with solving the normal equations directly. So &#8230; what do we do?<\/p>\n<p>Let&#8217;s break up the symmetric eigenvalue problem into these steps. With an input matrix <span class=\"math inline\">\\(A\\)<\/span>:<\/p>\n<ol style=\"list-style-type: decimal\">\n<li>Reduce <span class=\"math inline\">\\(A\\)<\/span> to tridiagonal <span class=\"math inline\">\\(T\\)<\/span>,<\/li>\n<li>Perform the first step of the <span class=\"math inline\">\\(QR\\)<\/span> iteration to find and apply <span class=\"math inline\">\\(Q_{q_1}\\)<\/span>,<\/li>\n<li>Perform the remainder of <span class=\"math inline\">\\(QR\\)<\/span> iteration until convergence, finding <span class=\"math inline\">\\(Q_{q_i}\\)<\/span> as in step 2.<\/li>\n<\/ol>\n<p>The total transformations on <span class=\"math inline\">\\(A\\)<\/span> come from reduction in step 1 and then repeating step 2 over and over. Let <span class=\"math inline\">\\(Q = \\prod Q_{q_i}\\ \\ Q_{Tri}\\)<\/span>. Then, <span class=\"math inline\">\\(Q^T A Q=D\\)<\/span>.<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>Here&#8217;s some equivalent steps for the SVD on an input matrix <span class=\"math inline\">\\(A\\)<\/span> (we need <span class=\"math inline\">\\(\\mathrm{eig}(A^TA)\\)<\/span>):<\/p>\n<ol style=\"list-style-type: decimal\">\n<li>Reduce <span class=\"math inline\">\\(A\\)<\/span> to <em>bidiagonal<\/em> <span class=\"math inline\">\\(B\\)<\/span> (note, <span class=\"math inline\">\\(B^TB\\)<\/span> is tridiagonal) using <span class=\"math inline\">\\(Q_{Bidiag}\\)<\/span>,<\/li>\n<li>Do some fancy footwork to compute <span class=\"math inline\">\\(Q_{q_1}\\)<\/span> as above (we only need parts of <span class=\"math inline\">\\(B^TB\\)<\/span>, not the whole thing) and then apply <span class=\"math inline\">\\(Q_{q_1}\\)<\/span> to <span class=\"math inline\">\\(B\\)<\/span>,<\/li>\n<li><span class=\"math inline\">\\(B Q_{q_1}\\)<\/span> is no longer bidiagonal. Use <code>bounce_blemish_down_diagonal<\/code> <a href=\"http:\/\/drsfenner.org\/blog\/2016\/03\/givens-rotations-and-the-case-of-the-blemished-bidiagonal-matrix\/\">(from the Givens\/Blemish post)<\/a> to restore it to blissful bidiagonality. Store or accumulate the products of the left <span class=\"math inline\">\\(G_{L_i}\\)<\/span> and right <span class=\"math inline\">\\(G_{R_i}\\)<\/span> Givens transformations that bouncing uses. You can see in <code>bounce_blemish_down_diagonal<\/code> that these come in <em>(left, right)<\/em> pairs followed by a last <em>left<\/em>,<\/li>\n<li>Perform the remainder of this iteration until convergence, finding <span class=\"math inline\">\\(Q_{q_i}\\)<\/span>, <span class=\"math inline\">\\(G_{L_i}\\)<\/span>, and <span class=\"math inline\">\\(G_{R_i}\\)<\/span> as in steps 2 and 3.<\/li>\n<\/ol>\n<p>The total transformations on <span class=\"math inline\">\\(B\\)<\/span> come from repeating steps 2 and 3 over and over. Let <span class=\"math inline\">\\(L = \\prod G_{L_i}\\)<\/span> and <span class=\"math inline\">\\(R = Q_{q_i} \\prod G_{R_i}\\)<\/span>. Then, more or less, <span class=\"math inline\">\\(L^T B R = D\\)<\/span>. Depending on your specific accounting\/accumulating and the return values you need, you might want different transposes in there. We often care about the form <span class=\"math inline\">\\(A=L_{all} D R^T_{all}\\)<\/span>. These are the values returned by, for example, LAPACK&#8217;s SVD routines (and hence, NumPy, SciPy, matlab, etc.). We get that by wrapping up (multiplying) the transformations that gave us <span class=\"math inline\">\\(A \\rightarrow B\\)<\/span> <em>and<\/em> the <span class=\"math inline\">\\(QR\\)<\/span>-step iterations into <span class=\"math inline\">\\(L_{all}\\)<\/span> and <span class=\"math inline\">\\(R_{all}\\)<\/span>. For example, <span class=\"math inline\">\\(L_{all}=L Q_{Bidiag}\\)<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3 id=\"by-hook-or-crook-essentially-similar\">By Hook or Crook: Essentially Similar<\/h3>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>By a wonderful theorem (the <em>Implicit Q Theorem<\/em>), since we started out with <span class=\"math inline\">\\(Q_{q_1}\\)<\/span> in both cases (<em>and<\/em> <span class=\"math inline\">\\(B\\)<\/span> and <span class=\"math inline\">\\(T\\)<\/span> were &quot;full&quot;\/unreduced <em>and<\/em> both processes end up with a tridiagonal &#8211; see G&amp;VL for details), the entire sequence of transformations (generated using <span class=\"math inline\">\\(A^TA\\)<\/span> in the QR process and <span class=\"math inline\">\\(A\\)<\/span>+fancy-footwork in the SVD process) will be &quot;essentially similar&quot; meaning they differ only by some signs on the diagonals.<\/p>\n<p>If we take <span class=\"math inline\">\\(T=B^TB\\)<\/span> at the start, then an update to <span class=\"math inline\">\\(T_{i+1}=Q^T_i T_i Q_i\\)<\/span> is related to <span class=\"math inline\">\\(B_i\\)<\/span> by applying the same <span class=\"math inline\">\\(Q_i\\)<\/span> to <span class=\"math inline\">\\(B_i\\)<\/span> and then &quot;fixing up&quot; <span class=\"math inline\">\\(B_i\\)<\/span> with an orthogonal <span class=\"math inline\">\\(U_i\\)<\/span> to maintain its bidiagonal form. So, <span class=\"math inline\">\\(B_{i+1}=U^T_i B_i Q_i\\)<\/span>.<\/p>\n<p><span class=\"math display\">\\[T_{i+1}=Q^T_i T_i Q_i = Q^T_i B^T_i B_i Q_i = Q^T_i B^T_i {U^T_i}^T U_i^T B_i Q_i<br \/>\n         =Q^T_i B^T_i U_i U^T_i B_i Q_i = B^T_{i+1} B_{i+1}\\]<\/span><\/p>\n<p>It&#8217;s pretty clear that that pairs <span class=\"math inline\">\\(L,R\\)<\/span> (by construction\/algorithm above) and <span class=\"math inline\">\\(Q,U\\)<\/span> (by &quot;proof of possibility&quot; &#8211; or an existence claim) are not directly related on a step-by-step basis. The existence claim is demonstrated by a series of equalities and it simply says <em>some<\/em> <span class=\"math inline\">\\(U,Q\\)<\/span> can hold that equality between <span class=\"math inline\">\\(B^T_iB_i\\)<\/span> and <span class=\"math inline\">\\(T_i\\)<\/span>. The ones we find with our QR and SVD iterations are different. Most importantly, our &quot;fixing up&quot; is not all in <span class=\"math inline\">\\(U\\)<\/span> (the orthogonal matrix on the interior of the expression); fixing up happens on <em>both<\/em> transformation, <span class=\"math inline\">\\(L\\)<\/span> and <span class=\"math inline\">\\(R\\)<\/span>. Regardless, once we gather them all up, they (<span class=\"math inline\">\\(L\\)<\/span> and <span class=\"math inline\">\\(Q\\)<\/span>) become essentially similar <em>and<\/em>, importantly, they lead to the correct <span class=\"math inline\">\\(D\\)<\/span>.<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3 id=\"roll-call\">Roll Call<\/h3>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>We are going to make use of many of the friends that we developed of the past several posts. For simplicity, I placed them (and an occasional modification or two) in <code>ghbs_utils.py<\/code>. That file is available at: http:\/\/drsfenner.org\/public\/notebooks\/additional\/ghbs_utils.py. Because of this import, I don&#8217;t believe this notebook will run (unmodified) on nbviewer.ipython.org. But, you can still download it and run it, see the link at the bottom of this post.<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[1]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"kn\">import<\/span> <span class=\"nn\">numpy<\/span> <span class=\"kn\">as<\/span> <span class=\"nn\">np<\/span>\n<span class=\"kn\">import<\/span> <span class=\"nn\">numpy.linalg<\/span> <span class=\"kn\">as<\/span> <span class=\"nn\">nla<\/span>\n<span class=\"kn\">from<\/span> <span class=\"nn\">pprint<\/span> <span class=\"kn\">import<\/span> <span class=\"n\">pprint<\/span>\n\n<span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">set_printoptions<\/span><span class=\"p\">(<\/span><span class=\"n\">suppress<\/span><span class=\"o\">=<\/span><span class=\"bp\">True<\/span><span class=\"p\">)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[2]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"c1\"># ghbs_utils.py available at:  <\/span>\n<span class=\"c1\"># http:\/\/drsfenner.org\/public\/notebooks\/additional\/ghbs_utils.py<\/span>\n<span class=\"kn\">from<\/span> <span class=\"nn\">ghbs_utils<\/span> <span class=\"kn\">import<\/span> <span class=\"p\">(<\/span><span class=\"n\">bounce_blemish_down_diagonal<\/span><span class=\"p\">,<\/span> <span class=\"n\">decouple_or_reduce<\/span><span class=\"p\">,<\/span>\n                        <span class=\"n\">zeroing_givens_coeffs<\/span><span class=\"p\">,<\/span> <span class=\"n\">left_givensT<\/span><span class=\"p\">,<\/span> <span class=\"n\">right_givens<\/span><span class=\"p\">,<\/span> \n                        <span class=\"n\">apply_givens_lefts_on_right<\/span><span class=\"p\">,<\/span> <span class=\"n\">apply_givens_rights_on_left<\/span><span class=\"p\">,<\/span>\n                        <span class=\"n\">make_upper_bidiag<\/span><span class=\"p\">,<\/span> <span class=\"n\">extract_upper_bidiag<\/span><span class=\"p\">,<\/span>\n                        <span class=\"n\">house_bidiag<\/span><span class=\"p\">,<\/span> <span class=\"n\">extract_packed_house_bidiag<\/span><span class=\"p\">)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3 id=\"svd-helpers-i-finding-last-non-zero-runs-in-a-sequence\">SVD Helpers I: Finding Last Non-Zero Runs in a Sequence<\/h3>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>When we attempt to apply the <span class=\"math inline\">\\(QR\\)<\/span>-iteration above, we have to have a &quot;full&quot; super-diagonal. We can find such a diagonal by looking for runs of non-zeros from the end of the super-diagonal backwards. The following code does that.<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[3]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"k\">def<\/span> <span class=\"nf\">find_first_run_boolean<\/span><span class=\"p\">(<\/span><span class=\"n\">vec<\/span><span class=\"p\">,<\/span> <span class=\"n\">front<\/span><span class=\"o\">=<\/span><span class=\"bp\">True<\/span><span class=\"p\">):<\/span>\n    <span class=\"sd\">&#39;&#39;&#39; use front=False if you will want a run from the end&#39;&#39;&#39;<\/span>\n    <span class=\"n\">lbl<\/span> <span class=\"o\">=<\/span> <span class=\"nb\">xrange<\/span><span class=\"p\">(<\/span><span class=\"nb\">len<\/span><span class=\"p\">(<\/span><span class=\"n\">vec<\/span><span class=\"p\">))<\/span>\n    <span class=\"k\">if<\/span> <span class=\"ow\">not<\/span> <span class=\"n\">front<\/span><span class=\"p\">:<\/span>\n        <span class=\"n\">vec<\/span> <span class=\"o\">=<\/span> <span class=\"nb\">reversed<\/span><span class=\"p\">(<\/span><span class=\"n\">vec<\/span><span class=\"p\">)<\/span>\n        <span class=\"n\">lbl<\/span> <span class=\"o\">=<\/span> <span class=\"nb\">reversed<\/span><span class=\"p\">(<\/span><span class=\"n\">lbl<\/span><span class=\"p\">)<\/span>\n    \n    <span class=\"n\">foundStart<\/span> <span class=\"o\">=<\/span> <span class=\"bp\">False<\/span>\n    <span class=\"k\">for<\/span> <span class=\"n\">idx<\/span><span class=\"p\">,<\/span> <span class=\"n\">itm<\/span> <span class=\"ow\">in<\/span> <span class=\"nb\">zip<\/span><span class=\"p\">(<\/span><span class=\"n\">lbl<\/span><span class=\"p\">,<\/span> <span class=\"n\">vec<\/span><span class=\"p\">):<\/span>\n        <span class=\"k\">if<\/span> <span class=\"n\">foundStart<\/span><span class=\"p\">:<\/span>\n            <span class=\"k\">if<\/span> <span class=\"n\">itm<\/span><span class=\"p\">:<\/span>\n                <span class=\"n\">end<\/span> <span class=\"o\">=<\/span> <span class=\"n\">idx<\/span>\n            <span class=\"k\">else<\/span><span class=\"p\">:<\/span>\n                <span class=\"k\">break<\/span>\n        <span class=\"k\">else<\/span><span class=\"p\">:<\/span>\n            <span class=\"k\">if<\/span> <span class=\"n\">itm<\/span><span class=\"p\">:<\/span>\n                <span class=\"n\">start<\/span> <span class=\"o\">=<\/span> <span class=\"n\">end<\/span> <span class=\"o\">=<\/span> <span class=\"n\">idx<\/span>\n                <span class=\"n\">foundStart<\/span> <span class=\"o\">=<\/span> <span class=\"bp\">True<\/span>\n    <span class=\"k\">if<\/span> <span class=\"ow\">not<\/span> <span class=\"n\">foundStart<\/span><span class=\"p\">:<\/span>\n        <span class=\"k\">return<\/span> <span class=\"nb\">slice<\/span><span class=\"p\">(<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">)<\/span> <span class=\"c1\"># equiv to slice(None, 0) [empty]<\/span>\n    <span class=\"k\">if<\/span> <span class=\"ow\">not<\/span> <span class=\"n\">front<\/span><span class=\"p\">:<\/span>\n        <span class=\"n\">start<\/span><span class=\"p\">,<\/span> <span class=\"n\">end<\/span> <span class=\"o\">=<\/span> <span class=\"n\">end<\/span><span class=\"p\">,<\/span> <span class=\"n\">start<\/span> <span class=\"c1\"># if we counted from end, swap<\/span>\n    <span class=\"k\">return<\/span> <span class=\"nb\">slice<\/span><span class=\"p\">(<\/span><span class=\"n\">start<\/span><span class=\"p\">,<\/span> <span class=\"n\">end<\/span><span class=\"o\">+<\/span><span class=\"mi\">1<\/span><span class=\"p\">)<\/span>  <span class=\"c1\"># allow for slicing use<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[4]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"n\">testcases<\/span> <span class=\"o\">=<\/span> <span class=\"p\">[(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">array<\/span><span class=\"p\">([<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span><span class=\"mi\">2<\/span><span class=\"p\">,<\/span><span class=\"mi\">3<\/span><span class=\"p\">,<\/span><span class=\"mi\">4<\/span><span class=\"p\">,<\/span><span class=\"mi\">5<\/span><span class=\"p\">,<\/span><span class=\"mi\">19<\/span><span class=\"p\">,<\/span><span class=\"mi\">20<\/span><span class=\"p\">]),<\/span> <span class=\"nb\">slice<\/span><span class=\"p\">(<\/span><span class=\"mi\">2<\/span><span class=\"p\">,<\/span><span class=\"mi\">9<\/span><span class=\"p\">,<\/span><span class=\"bp\">None<\/span><span class=\"p\">)),<\/span>\n             <span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">zeros<\/span><span class=\"p\">(<\/span><span class=\"mi\">5<\/span><span class=\"p\">),<\/span> <span class=\"nb\">slice<\/span><span class=\"p\">(<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"bp\">None<\/span><span class=\"p\">)),<\/span>\n             <span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">array<\/span><span class=\"p\">([<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span><span class=\"mi\">2<\/span><span class=\"p\">,<\/span><span class=\"mi\">3<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">4<\/span><span class=\"p\">,<\/span><span class=\"mi\">5<\/span><span class=\"p\">]),<\/span> <span class=\"nb\">slice<\/span><span class=\"p\">(<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">3<\/span><span class=\"p\">,<\/span><span class=\"bp\">None<\/span><span class=\"p\">)),<\/span>\n             <span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">array<\/span><span class=\"p\">([<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">10<\/span><span class=\"p\">]),<\/span> <span class=\"nb\">slice<\/span><span class=\"p\">(<\/span><span class=\"mi\">4<\/span><span class=\"p\">,<\/span><span class=\"mi\">5<\/span><span class=\"p\">,<\/span><span class=\"bp\">None<\/span><span class=\"p\">)),<\/span>\n            <span class=\"p\">]<\/span>\n\n<span class=\"k\">print<\/span> <span class=\"nb\">all<\/span><span class=\"p\">(<\/span><span class=\"n\">t_out<\/span> <span class=\"o\">==<\/span> <span class=\"n\">find_first_run_boolean<\/span><span class=\"p\">(<\/span><span class=\"n\">t_in<\/span><span class=\"p\">,<\/span> <span class=\"n\">front<\/span><span class=\"o\">=<\/span><span class=\"bp\">True<\/span><span class=\"p\">)<\/span> <span class=\"k\">for<\/span> <span class=\"n\">t_in<\/span><span class=\"p\">,<\/span> <span class=\"n\">t_out<\/span> <span class=\"ow\">in<\/span> <span class=\"n\">testcases<\/span><span class=\"p\">)<\/span>\n<span class=\"c1\">#for t_input, t_output in testcases:<\/span>\n<span class=\"c1\">#    sli = find_first_run_boolean(t_input, front=True)<\/span>\n<span class=\"c1\">#    print sli == t_output<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"output_wrapper\">\n<div class=\"output\">\n<div class=\"output_area\">\n<div class=\"prompt\"><\/div>\n<div class=\"output_subarea output_stream output_stdout output_text\">\n<pre>True\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[5]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"n\">testcases<\/span> <span class=\"o\">=<\/span> <span class=\"p\">[(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">array<\/span><span class=\"p\">([<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span><span class=\"mi\">2<\/span><span class=\"p\">,<\/span><span class=\"mi\">3<\/span><span class=\"p\">,<\/span><span class=\"mi\">4<\/span><span class=\"p\">,<\/span><span class=\"mi\">5<\/span><span class=\"p\">,<\/span><span class=\"mi\">19<\/span><span class=\"p\">,<\/span><span class=\"mi\">20<\/span><span class=\"p\">]),<\/span> <span class=\"nb\">slice<\/span><span class=\"p\">(<\/span><span class=\"mi\">2<\/span><span class=\"p\">,<\/span><span class=\"mi\">9<\/span><span class=\"p\">,<\/span><span class=\"bp\">None<\/span><span class=\"p\">)),<\/span>\n             <span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">zeros<\/span><span class=\"p\">(<\/span><span class=\"mi\">5<\/span><span class=\"p\">),<\/span> <span class=\"nb\">slice<\/span><span class=\"p\">(<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"bp\">None<\/span><span class=\"p\">)),<\/span>\n             <span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">array<\/span><span class=\"p\">([<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span><span class=\"mi\">2<\/span><span class=\"p\">,<\/span><span class=\"mi\">3<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">4<\/span><span class=\"p\">,<\/span><span class=\"mi\">5<\/span><span class=\"p\">]),<\/span> <span class=\"nb\">slice<\/span><span class=\"p\">(<\/span><span class=\"mi\">5<\/span><span class=\"p\">,<\/span><span class=\"mi\">7<\/span><span class=\"p\">,<\/span><span class=\"bp\">None<\/span><span class=\"p\">)),<\/span>\n             <span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">array<\/span><span class=\"p\">([<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">10<\/span><span class=\"p\">]),<\/span> <span class=\"nb\">slice<\/span><span class=\"p\">(<\/span><span class=\"mi\">4<\/span><span class=\"p\">,<\/span><span class=\"mi\">5<\/span><span class=\"p\">,<\/span><span class=\"bp\">None<\/span><span class=\"p\">)),<\/span>\n            <span class=\"p\">]<\/span>\n<span class=\"k\">print<\/span> <span class=\"nb\">all<\/span><span class=\"p\">(<\/span><span class=\"n\">t_out<\/span> <span class=\"o\">==<\/span> <span class=\"n\">find_first_run_boolean<\/span><span class=\"p\">(<\/span><span class=\"n\">t_in<\/span><span class=\"p\">,<\/span> <span class=\"n\">front<\/span><span class=\"o\">=<\/span><span class=\"bp\">False<\/span><span class=\"p\">)<\/span> <span class=\"k\">for<\/span> <span class=\"n\">t_in<\/span><span class=\"p\">,<\/span> <span class=\"n\">t_out<\/span> <span class=\"ow\">in<\/span> <span class=\"n\">testcases<\/span><span class=\"p\">)<\/span>\n<span class=\"c1\">#for t_input, t_output in testcases:<\/span>\n<span class=\"c1\">#    sli = find_first_run_boolean(t_input, front=False)<\/span>\n<span class=\"c1\">#    print sli == t_output<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"output_wrapper\">\n<div class=\"output\">\n<div class=\"output_area\">\n<div class=\"prompt\"><\/div>\n<div class=\"output_subarea output_stream output_stdout output_text\">\n<pre>True\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3 id=\"svd-helpers-ii-trimming-tiny-values\">SVD Helpers II: Trimming Tiny Values<\/h3>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>There are several possibilities for determining when a value is numerically small enough to consider it &quot;practically&quot; zero. G&amp;VL do so by comparing an super-diagonal element to its diagonal &quot;friend&quot; (on diagonal in the same row) and the diagonal element in the row below it. I also use this opportunity to remove small values on the diagonal itself (even though G&amp;VL don&#8217;t mention that step in their algorithm).<\/p>\n<p>One trick that this code uses is the <code>ndarray.flat<\/code> iterator. This is a view into the entire array as a flat entity. By starting at index 1 and taking steps of in the size of the number of elements in a row plus one, we get each super-diagonal entry. There is a NumPy function to do similar for extraction (<code>np.diag<\/code>), but it can&#8217;t be used as a target for an assignment. So, I use <code>ndarray.flat<\/code> several times to give a feel for what it is doing. A side note: there is also <code>np.diag_indices<\/code> which I really dislike, but your mileage might vary!<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[6]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"c1\"># this is a factored-out function for the <\/span>\n<span class=\"c1\"># first loop step in G&amp;VL Algo. 8.6.2<\/span>\n<span class=\"k\">def<\/span> <span class=\"nf\">remove_small<\/span><span class=\"p\">(<\/span><span class=\"n\">B<\/span><span class=\"p\">,<\/span> <span class=\"n\">eps<\/span><span class=\"p\">):<\/span>\n    <span class=\"n\">m<\/span><span class=\"p\">,<\/span><span class=\"n\">n<\/span> <span class=\"o\">=<\/span> <span class=\"n\">B<\/span><span class=\"o\">.<\/span><span class=\"n\">shape<\/span>\n    <span class=\"n\">diag<\/span>          <span class=\"o\">=<\/span> <span class=\"n\">B<\/span><span class=\"o\">.<\/span><span class=\"n\">flat<\/span><span class=\"p\">[<\/span> <span class=\"p\">::<\/span><span class=\"n\">n<\/span><span class=\"o\">+<\/span><span class=\"mi\">1<\/span><span class=\"p\">]<\/span> <span class=\"c1\"># np.diag is fine, but want to explain B.flat below<\/span>\n    <span class=\"n\">abs_diag<\/span>      <span class=\"o\">=<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">abs<\/span><span class=\"p\">(<\/span><span class=\"n\">diag<\/span><span class=\"p\">)<\/span>\n    <span class=\"n\">B<\/span><span class=\"o\">.<\/span><span class=\"n\">flat<\/span><span class=\"p\">[::<\/span><span class=\"n\">n<\/span><span class=\"o\">+<\/span><span class=\"mi\">1<\/span><span class=\"p\">]<\/span> <span class=\"o\">=<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">where<\/span><span class=\"p\">(<\/span><span class=\"n\">abs_diag<\/span> <span class=\"o\">&lt;<\/span> <span class=\"n\">eps<\/span><span class=\"p\">,<\/span> <span class=\"mf\">0.0<\/span><span class=\"p\">,<\/span> <span class=\"n\">diag<\/span><span class=\"p\">)<\/span>\n\n    <span class=\"n\">super_diag<\/span>     <span class=\"o\">=<\/span> <span class=\"n\">B<\/span><span class=\"o\">.<\/span><span class=\"n\">flat<\/span><span class=\"p\">[<\/span><span class=\"mi\">1<\/span><span class=\"p\">::<\/span><span class=\"n\">n<\/span><span class=\"o\">+<\/span><span class=\"mi\">1<\/span><span class=\"p\">]<\/span>\n    \n    <span class=\"c1\"># superdiag -&gt; 0 if superdiag &lt;= eps * (this row diag + next row diag)<\/span>\n    <span class=\"c1\"># where eps is small multiple of unit roundoff<\/span>\n    <span class=\"n\">small_val<\/span> <span class=\"o\">=<\/span> <span class=\"n\">eps<\/span> <span class=\"o\">*<\/span> <span class=\"p\">(<\/span><span class=\"n\">abs_diag<\/span><span class=\"p\">[:<\/span><span class=\"o\">-<\/span><span class=\"mi\">1<\/span><span class=\"p\">]<\/span> <span class=\"o\">+<\/span> <span class=\"n\">abs_diag<\/span><span class=\"p\">[<\/span><span class=\"mi\">1<\/span><span class=\"p\">:])<\/span> <span class=\"c1\"># or just = eps<\/span>\n    <span class=\"n\">B<\/span><span class=\"o\">.<\/span><span class=\"n\">flat<\/span><span class=\"p\">[<\/span><span class=\"mi\">1<\/span><span class=\"p\">::<\/span><span class=\"n\">n<\/span><span class=\"o\">+<\/span><span class=\"mi\">1<\/span><span class=\"p\">]<\/span> <span class=\"o\">=<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">where<\/span><span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">abs<\/span><span class=\"p\">(<\/span><span class=\"n\">super_diag<\/span><span class=\"p\">)<\/span> <span class=\"o\">&lt;<\/span> <span class=\"n\">small_val<\/span><span class=\"p\">,<\/span> <span class=\"mf\">0.0<\/span><span class=\"p\">,<\/span> <span class=\"n\">super_diag<\/span><span class=\"p\">)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3 id=\"svd-helpers-iii-cleaning-and-partitioning\">SVD Helpers III: Cleaning and Partitioning<\/h3>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>Here we have the code that is responsible for find the biggest unsolved subproblem towards the lower-right corner of our input matrix. Everything further to the lower-right is already reduced to a pure diagonal; hence it already has its singular values computed. Everything to the upper-left still needs to be examined.<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[7]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"c1\"># this is a factored-out function for the <\/span>\n<span class=\"c1\"># first two steps in G&amp;VL Algo. 8.6.2<\/span>\n<span class=\"k\">def<\/span> <span class=\"nf\">clean_and_partition<\/span><span class=\"p\">(<\/span><span class=\"n\">B<\/span><span class=\"p\">,<\/span> <span class=\"n\">eps<\/span><span class=\"o\">=.<\/span><span class=\"mo\">00001<\/span><span class=\"p\">):<\/span>\n    <span class=\"sd\">&#39;&#39;&#39;return both the sub-matrix and the index it starts at&#39;&#39;&#39;<\/span>\n    <span class=\"n\">remove_small<\/span><span class=\"p\">(<\/span><span class=\"n\">B<\/span><span class=\"p\">,<\/span> <span class=\"n\">eps<\/span><span class=\"o\">=<\/span><span class=\"n\">eps<\/span><span class=\"p\">)<\/span>\n    \n    <span class=\"n\">m<\/span><span class=\"p\">,<\/span><span class=\"n\">n<\/span> <span class=\"o\">=<\/span> <span class=\"n\">B<\/span><span class=\"o\">.<\/span><span class=\"n\">shape<\/span>\n    <span class=\"n\">super_diag<\/span> <span class=\"o\">=<\/span> <span class=\"n\">B<\/span><span class=\"o\">.<\/span><span class=\"n\">flat<\/span><span class=\"p\">[<\/span><span class=\"mi\">1<\/span><span class=\"p\">::<\/span><span class=\"n\">n<\/span><span class=\"o\">+<\/span><span class=\"mi\">1<\/span><span class=\"p\">]<\/span>\n    <span class=\"n\">super_slice<\/span> <span class=\"o\">=<\/span> <span class=\"n\">find_first_run_boolean<\/span><span class=\"p\">(<\/span><span class=\"n\">super_diag<\/span><span class=\"p\">,<\/span> <span class=\"n\">front<\/span><span class=\"o\">=<\/span><span class=\"bp\">False<\/span><span class=\"p\">)<\/span>\n    <span class=\"k\">if<\/span> <span class=\"ow\">not<\/span> <span class=\"n\">super_slice<\/span><span class=\"o\">.<\/span><span class=\"n\">stop<\/span><span class=\"p\">:<\/span>\n        <span class=\"k\">return<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">empty<\/span><span class=\"p\">(<\/span><span class=\"mi\">0<\/span><span class=\"p\">),<\/span> <span class=\"mi\">0<\/span>\n        \n    <span class=\"c1\"># based on the super diagonal, we want <\/span>\n    <span class=\"n\">full_slice<\/span> <span class=\"o\">=<\/span> <span class=\"nb\">slice<\/span><span class=\"p\">(<\/span><span class=\"n\">super_slice<\/span><span class=\"o\">.<\/span><span class=\"n\">start<\/span><span class=\"p\">,<\/span> <span class=\"n\">super_slice<\/span><span class=\"o\">.<\/span><span class=\"n\">stop<\/span> <span class=\"o\">+<\/span> <span class=\"mi\">1<\/span><span class=\"p\">)<\/span>\n    <span class=\"k\">return<\/span> <span class=\"n\">B<\/span><span class=\"p\">[<\/span><span class=\"n\">full_slice<\/span><span class=\"p\">,<\/span> <span class=\"n\">full_slice<\/span><span class=\"p\">],<\/span> <span class=\"n\">super_slice<\/span><span class=\"o\">.<\/span><span class=\"n\">start<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[8]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"n\">ubd<\/span> <span class=\"o\">=<\/span> <span class=\"n\">make_upper_bidiag<\/span><span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">arange<\/span><span class=\"p\">(<\/span><span class=\"mi\">8<\/span><span class=\"p\">)<\/span><span class=\"o\">*<\/span><span class=\"mi\">2<\/span> <span class=\"o\">+<\/span> <span class=\"mi\">10<\/span><span class=\"p\">,<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">arange<\/span><span class=\"p\">(<\/span><span class=\"mi\">7<\/span><span class=\"p\">)<\/span><span class=\"o\">*<\/span><span class=\"mi\">2<\/span><span class=\"o\">+<\/span><span class=\"mi\">1<\/span><span class=\"p\">)<\/span>\n<span class=\"n\">ubd<\/span><span class=\"p\">[<\/span><span class=\"mi\">3<\/span><span class=\"p\">,<\/span><span class=\"mi\">4<\/span><span class=\"p\">]<\/span> <span class=\"o\">=<\/span> <span class=\"mi\">0<\/span>\n<span class=\"k\">print<\/span> <span class=\"n\">ubd<\/span>\n\n<span class=\"n\">B<\/span><span class=\"p\">,<\/span><span class=\"n\">k<\/span> <span class=\"o\">=<\/span> <span class=\"n\">clean_and_partition<\/span><span class=\"p\">(<\/span><span class=\"n\">ubd<\/span><span class=\"p\">,<\/span> <span class=\"o\">.<\/span><span class=\"mo\">01<\/span><span class=\"p\">)<\/span>\n<span class=\"k\">print<\/span> <span class=\"n\">B<\/span>\n<span class=\"k\">print<\/span> <span class=\"n\">k<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"output_wrapper\">\n<div class=\"output\">\n<div class=\"output_area\">\n<div class=\"prompt\"><\/div>\n<div class=\"output_subarea output_stream output_stdout output_text\">\n<pre>[[10  1  0  0  0  0  0  0]\n [ 0 12  3  0  0  0  0  0]\n [ 0  0 14  5  0  0  0  0]\n [ 0  0  0 16  0  0  0  0]\n [ 0  0  0  0 18  9  0  0]\n [ 0  0  0  0  0 20 11  0]\n [ 0  0  0  0  0  0 22 13]\n [ 0  0  0  0  0  0  0 24]]\n[[18  9  0  0]\n [ 0 20 11  0]\n [ 0  0 22 13]\n [ 0  0  0 24]]\n4\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3 id=\"pushing-towards-diagonality\">Pushing Towards Diagonality<\/h3>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>One thing I completely avoided discussing above is the notion of the &quot;Wilkinson shift&quot;. Essentially, this is a technique that speeds the convergence of the <span class=\"math inline\">\\(QR\\)<\/span>-iteration. See G&amp;VL, 3rd, pg. 418-420 for details. All it does is find an eigenvalue of the lower-right <span class=\"math inline\">\\(2\\ x\\ 2\\)<\/span> matrix that is closer to the bottom-right element and subtract that from our diagonal entries.<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[9]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"c1\"># G&amp;VL algorithm 8.6.1 on page ~ 454 .. <\/span>\n<span class=\"k\">def<\/span> <span class=\"nf\">step_full_bidiagonal_towards_diagonal<\/span><span class=\"p\">(<\/span><span class=\"n\">B<\/span><span class=\"p\">):<\/span>\n    <span class=\"c1\"># let T &lt;- trailing 2x2 submatrix of B.T.dot(B)<\/span>\n    <span class=\"c1\"># T = B.T.dot(B)[-2:,-2:]  .... row by col rule means we<\/span>\n    <span class=\"c1\">#     need last two rows of B.T (which are last two cols of B)<\/span>\n    <span class=\"c1\">#     and last two cols of B (which means we just need ....)<\/span>\n    <span class=\"n\">lastCols<\/span> <span class=\"o\">=<\/span> <span class=\"n\">B<\/span><span class=\"p\">[:,<\/span> <span class=\"o\">-<\/span><span class=\"mi\">2<\/span><span class=\"p\">:]<\/span>\n    <span class=\"n\">T<\/span> <span class=\"o\">=<\/span> <span class=\"n\">lastCols<\/span><span class=\"o\">.<\/span><span class=\"n\">T<\/span><span class=\"o\">.<\/span><span class=\"n\">dot<\/span><span class=\"p\">(<\/span><span class=\"n\">lastCols<\/span><span class=\"p\">)<\/span>\n    \n    <span class=\"c1\"># compute wilkinson shift value (&quot;mu&quot;)<\/span>\n    <span class=\"n\">d<\/span>  <span class=\"o\">=<\/span> <span class=\"p\">(<\/span><span class=\"n\">T<\/span><span class=\"p\">[<\/span><span class=\"o\">-<\/span><span class=\"mi\">2<\/span><span class=\"p\">,<\/span><span class=\"o\">-<\/span><span class=\"mi\">2<\/span><span class=\"p\">]<\/span> <span class=\"o\">-<\/span> <span class=\"n\">T<\/span><span class=\"p\">[<\/span><span class=\"o\">-<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span><span class=\"o\">-<\/span><span class=\"mi\">1<\/span><span class=\"p\">])<\/span> <span class=\"o\">\/<\/span> <span class=\"mf\">2.0<\/span>\n    <span class=\"n\">shift<\/span> <span class=\"o\">=<\/span> <span class=\"n\">T<\/span><span class=\"p\">[<\/span><span class=\"o\">-<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span><span class=\"o\">-<\/span><span class=\"mi\">1<\/span><span class=\"p\">]<\/span> <span class=\"o\">-<\/span> <span class=\"n\">T<\/span><span class=\"p\">[<\/span><span class=\"o\">-<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span><span class=\"o\">-<\/span><span class=\"mi\">2<\/span><span class=\"p\">]<\/span><span class=\"o\">**<\/span><span class=\"mi\">2<\/span> <span class=\"o\">\/<\/span> <span class=\"p\">(<\/span><span class=\"n\">d<\/span> <span class=\"o\">+<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">sign<\/span><span class=\"p\">(<\/span><span class=\"n\">d<\/span><span class=\"p\">)<\/span><span class=\"o\">*<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">sqrt<\/span><span class=\"p\">(<\/span><span class=\"n\">d<\/span><span class=\"o\">**<\/span><span class=\"mi\">2<\/span><span class=\"o\">+<\/span><span class=\"n\">T<\/span><span class=\"p\">[<\/span><span class=\"o\">-<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span><span class=\"o\">-<\/span><span class=\"mi\">2<\/span><span class=\"p\">]<\/span><span class=\"o\">**<\/span><span class=\"mi\">2<\/span><span class=\"p\">))<\/span>\n    \n    <span class=\"c1\"># special case givens_zero_right(B) on shift<\/span>\n    <span class=\"c1\"># the only explict part of T we need in this reduction step<\/span>\n    <span class=\"n\">T_00<\/span> <span class=\"o\">=<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">sum<\/span><span class=\"p\">(<\/span><span class=\"n\">B<\/span><span class=\"o\">.<\/span><span class=\"n\">T<\/span><span class=\"p\">[<\/span><span class=\"mi\">0<\/span><span class=\"p\">,:]<\/span><span class=\"o\">**<\/span><span class=\"mi\">2<\/span><span class=\"p\">)<\/span>\n    <span class=\"n\">T_01<\/span> <span class=\"o\">=<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">dot<\/span><span class=\"p\">(<\/span><span class=\"n\">B<\/span><span class=\"o\">.<\/span><span class=\"n\">T<\/span><span class=\"p\">[<\/span><span class=\"mi\">0<\/span><span class=\"p\">,:],<\/span> <span class=\"n\">B<\/span><span class=\"p\">[:,<\/span><span class=\"mi\">1<\/span><span class=\"p\">])<\/span> <span class=\"c1\">#B^T B --&gt; B.T[0,:], B[:,1] <\/span>\n\n    <span class=\"n\">c<\/span><span class=\"p\">,<\/span><span class=\"n\">s<\/span> <span class=\"o\">=<\/span> <span class=\"n\">zeroing_givens_coeffs<\/span><span class=\"p\">(<\/span><span class=\"n\">T_00<\/span><span class=\"o\">-<\/span><span class=\"n\">shift<\/span><span class=\"p\">,<\/span> <span class=\"n\">T_01<\/span><span class=\"p\">)<\/span>\n    <span class=\"n\">local_right<\/span> <span class=\"o\">=<\/span> <span class=\"p\">(<\/span><span class=\"n\">c<\/span><span class=\"p\">,<\/span><span class=\"n\">s<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">1<\/span><span class=\"p\">)<\/span> <span class=\"c1\"># store for later<\/span>\n    <span class=\"n\">right_givens<\/span><span class=\"p\">(<\/span><span class=\"n\">B<\/span><span class=\"p\">,<\/span> <span class=\"n\">local_right<\/span><span class=\"p\">)<\/span>\n    \n    <span class=\"n\">lefts<\/span><span class=\"p\">,<\/span> <span class=\"n\">rights<\/span> <span class=\"o\">=<\/span> <span class=\"n\">bounce_blemish_down_diagonal<\/span><span class=\"p\">(<\/span><span class=\"n\">B<\/span><span class=\"p\">)<\/span>\n    <span class=\"n\">rights<\/span> <span class=\"o\">=<\/span> <span class=\"p\">[<\/span><span class=\"n\">local_right<\/span><span class=\"p\">]<\/span> <span class=\"o\">+<\/span> <span class=\"n\">rights<\/span> <span class=\"c1\"># BONUS:  remove the copy (deque or reverse order)<\/span>\n    <span class=\"k\">return<\/span> <span class=\"n\">lefts<\/span><span class=\"p\">,<\/span> <span class=\"n\">rights<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>Note the order we apply <code>lefts<\/code> and <code>rights<\/code>. Since we are accumulating them &quot;outside-in&quot;, we need to apply them &quot;looking out&quot; to the accumulating <code>U<\/code> and <code>V<\/code>. If we applied them to <code>A<\/code>, we could do them in the expected order, but we wouldn&#8217;t have separate left and right transformations. <code>apply_givens_lefts_on_right<\/code> and its friend take care of these details. They are simply an accumulator loop starting with an identity matrix and multiplying in Givens rotations.<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[10]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"n\">A<\/span> <span class=\"o\">=<\/span> <span class=\"n\">make_upper_bidiag<\/span><span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">arange<\/span><span class=\"p\">(<\/span><span class=\"mf\">4.0<\/span><span class=\"p\">)<\/span><span class=\"o\">+<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">ones<\/span><span class=\"p\">(<\/span><span class=\"mi\">4<\/span><span class=\"p\">))<\/span>\n<span class=\"n\">A_orig<\/span> <span class=\"o\">=<\/span> <span class=\"n\">A<\/span><span class=\"o\">.<\/span><span class=\"n\">copy<\/span><span class=\"p\">()<\/span>\n\n<span class=\"n\">L<\/span><span class=\"p\">,<\/span>  <span class=\"n\">R<\/span>  <span class=\"o\">=<\/span> <span class=\"n\">step_full_bidiagonal_towards_diagonal<\/span><span class=\"p\">(<\/span><span class=\"n\">A<\/span><span class=\"p\">)<\/span>\n<span class=\"n\">L2<\/span><span class=\"p\">,<\/span> <span class=\"n\">R2<\/span> <span class=\"o\">=<\/span> <span class=\"n\">step_full_bidiagonal_towards_diagonal<\/span><span class=\"p\">(<\/span><span class=\"n\">A<\/span><span class=\"p\">)<\/span>\n\n<span class=\"c1\"># gather rotations into U and V<\/span>\n<span class=\"n\">L<\/span><span class=\"o\">.<\/span><span class=\"n\">extend<\/span><span class=\"p\">(<\/span><span class=\"n\">L2<\/span><span class=\"p\">)<\/span>\n<span class=\"n\">R<\/span><span class=\"o\">.<\/span><span class=\"n\">extend<\/span><span class=\"p\">(<\/span><span class=\"n\">R2<\/span><span class=\"p\">)<\/span>\n<span class=\"n\">U<\/span><span class=\"p\">,<\/span><span class=\"n\">V<\/span> <span class=\"o\">=<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">eye<\/span><span class=\"p\">(<\/span><span class=\"mi\">4<\/span><span class=\"p\">),<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">eye<\/span><span class=\"p\">(<\/span><span class=\"mi\">4<\/span><span class=\"p\">)<\/span>\n<span class=\"n\">apply_givens_lefts_on_right<\/span><span class=\"p\">(<\/span><span class=\"n\">U<\/span><span class=\"p\">,<\/span> <span class=\"n\">L<\/span><span class=\"p\">,<\/span> <span class=\"mi\">0<\/span><span class=\"p\">)<\/span>\n<span class=\"n\">apply_givens_rights_on_left<\/span><span class=\"p\">(<\/span><span class=\"n\">R<\/span><span class=\"p\">,<\/span> <span class=\"n\">V<\/span><span class=\"p\">,<\/span> <span class=\"mi\">0<\/span><span class=\"p\">)<\/span>\n\n<span class=\"c1\"># seems to be good enough<\/span>\n<span class=\"k\">print<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">allclose<\/span><span class=\"p\">(<\/span><span class=\"n\">A_orig<\/span><span class=\"p\">,<\/span> <span class=\"n\">U<\/span><span class=\"o\">.<\/span><span class=\"n\">dot<\/span><span class=\"p\">(<\/span><span class=\"n\">A<\/span><span class=\"p\">)<\/span><span class=\"o\">.<\/span><span class=\"n\">dot<\/span><span class=\"p\">(<\/span><span class=\"n\">V<\/span><span class=\"p\">))<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"output_wrapper\">\n<div class=\"output\">\n<div class=\"output_area\">\n<div class=\"prompt\"><\/div>\n<div class=\"output_subarea output_stream output_stdout output_text\">\n<pre>True\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3 id=\"and-full-svd\">And Full SVD<\/h3>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>With all the pieces in order, the find SVD algorithm isn&#8217;t too long. Having names for the sub-pieces, also gives us a mental break trying to interpret what is happening at each step. I left a few commented <code>assert<\/code>s in the code. It is a nice mental reminder that we are maintaining the equality wiht the original input throughout this process. Salud!<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[11]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"c1\"># GvL    page 454 - 456  algo 8.6.2 (full svd)<\/span>\n<span class=\"k\">def<\/span> <span class=\"nf\">clear_svd<\/span><span class=\"p\">(<\/span><span class=\"n\">A<\/span><span class=\"p\">,<\/span> <span class=\"n\">eps<\/span><span class=\"o\">=<\/span><span class=\"mf\">0.001<\/span><span class=\"p\">):<\/span> <span class=\"c1\"># eps=np.finfo(np.double).eps):<\/span>\n    <span class=\"n\">m<\/span><span class=\"p\">,<\/span><span class=\"n\">n<\/span> <span class=\"o\">=<\/span> <span class=\"n\">A<\/span><span class=\"o\">.<\/span><span class=\"n\">shape<\/span>\n    \n    <span class=\"c1\"># for safety check:  A_orig = A.copy()<\/span>\n    <span class=\"n\">both_betas<\/span> <span class=\"o\">=<\/span> <span class=\"n\">house_bidiag<\/span><span class=\"p\">(<\/span><span class=\"n\">A<\/span><span class=\"p\">)<\/span> <span class=\"c1\"># betas for U and Vt<\/span>\n    \n    <span class=\"c1\"># this sends back V; not worrying about it<\/span>\n    <span class=\"n\">U<\/span><span class=\"p\">,<\/span> <span class=\"n\">B<\/span><span class=\"p\">,<\/span> <span class=\"n\">V<\/span> <span class=\"o\">=<\/span> <span class=\"n\">extract_packed_house_bidiag<\/span><span class=\"p\">(<\/span><span class=\"n\">A<\/span><span class=\"p\">,<\/span> <span class=\"o\">*<\/span><span class=\"n\">both_betas<\/span><span class=\"p\">)<\/span>\n    <span class=\"c1\"># check:  assert np.allclose(A_orig, U.dot(B).dot(V.T))    <\/span>\n    <span class=\"n\">Vt<\/span> <span class=\"o\">=<\/span> <span class=\"n\">V<\/span><span class=\"o\">.<\/span><span class=\"n\">T<\/span> \n    \n    <span class=\"c1\"># B22 is the bottom-right-most &quot;full bidiagonal&quot; square block matrix of B<\/span>\n    <span class=\"c1\"># it is a view into B, so when B22 is updated, so is B<\/span>\n    <span class=\"n\">B22<\/span><span class=\"p\">,<\/span> <span class=\"n\">k<\/span> <span class=\"o\">=<\/span> <span class=\"n\">clean_and_partition<\/span><span class=\"p\">(<\/span><span class=\"n\">B<\/span><span class=\"p\">,<\/span><span class=\"n\">eps<\/span><span class=\"o\">=<\/span><span class=\"n\">eps<\/span><span class=\"p\">)<\/span>\n    <span class=\"k\">while<\/span> <span class=\"n\">B22<\/span><span class=\"o\">.<\/span><span class=\"n\">size<\/span><span class=\"p\">:<\/span>\n        <span class=\"n\">zero_idxs<\/span> <span class=\"o\">=<\/span> <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">where<\/span><span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">abs<\/span><span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">diag<\/span><span class=\"p\">(<\/span><span class=\"n\">B22<\/span><span class=\"p\">))<\/span> <span class=\"o\">&lt;<\/span> <span class=\"n\">eps<\/span><span class=\"p\">)[<\/span><span class=\"mi\">0<\/span><span class=\"p\">]<\/span>\n        <span class=\"k\">if<\/span> <span class=\"n\">zero_idxs<\/span><span class=\"o\">.<\/span><span class=\"n\">size<\/span><span class=\"p\">:<\/span>\n            <span class=\"n\">lefts<\/span><span class=\"p\">,<\/span> <span class=\"n\">rights<\/span> <span class=\"o\">=<\/span> <span class=\"n\">decouple_or_reduce<\/span><span class=\"p\">(<\/span><span class=\"n\">B22<\/span><span class=\"p\">,<\/span> <span class=\"n\">zero_idxs<\/span><span class=\"p\">)<\/span>\n        <span class=\"k\">else<\/span><span class=\"p\">:<\/span>\n            <span class=\"n\">lefts<\/span><span class=\"p\">,<\/span> <span class=\"n\">rights<\/span> <span class=\"o\">=<\/span> <span class=\"n\">step_full_bidiagonal_towards_diagonal<\/span><span class=\"p\">(<\/span><span class=\"n\">B22<\/span><span class=\"p\">)<\/span>\n        \n        <span class=\"c1\"># this feels weird, but since U and Vt are accumulating the effects on<\/span>\n        <span class=\"c1\"># B22, we develop them on one side (on B22) and apply them on the respective<\/span>\n        <span class=\"c1\"># (outer) transformation (U or Vt)<\/span>\n        <span class=\"n\">apply_givens_lefts_on_right<\/span><span class=\"p\">(<\/span><span class=\"n\">U<\/span><span class=\"p\">,<\/span> <span class=\"n\">lefts<\/span><span class=\"p\">,<\/span>   <span class=\"n\">shift_dex<\/span><span class=\"o\">=<\/span><span class=\"n\">k<\/span><span class=\"p\">)<\/span>\n        <span class=\"n\">apply_givens_rights_on_left<\/span><span class=\"p\">(<\/span><span class=\"n\">rights<\/span><span class=\"p\">,<\/span> <span class=\"n\">Vt<\/span><span class=\"p\">,<\/span> <span class=\"n\">shift_dex<\/span><span class=\"o\">=<\/span><span class=\"n\">k<\/span><span class=\"p\">)<\/span>\n\n        <span class=\"c1\"># check:  assert np.allclose(A_orig, U.dot(B).dot(Vt))<\/span>\n        \n        <span class=\"n\">B22<\/span><span class=\"p\">,<\/span> <span class=\"n\">k<\/span> <span class=\"o\">=<\/span> <span class=\"n\">clean_and_partition<\/span><span class=\"p\">(<\/span><span class=\"n\">B<\/span><span class=\"p\">,<\/span><span class=\"n\">eps<\/span><span class=\"o\">=<\/span><span class=\"n\">eps<\/span><span class=\"p\">)<\/span>\n    <span class=\"k\">return<\/span> <span class=\"n\">U<\/span><span class=\"p\">,<\/span> <span class=\"n\">B<\/span><span class=\"p\">,<\/span> <span class=\"n\">Vt<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3 id=\"some-tests\">Some Tests<\/h3>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[12]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"k\">def<\/span> <span class=\"nf\">test_check<\/span><span class=\"p\">(<\/span><span class=\"n\">M<\/span><span class=\"p\">):<\/span>\n    <span class=\"n\">expected_S<\/span> <span class=\"o\">=<\/span> <span class=\"n\">nla<\/span><span class=\"o\">.<\/span><span class=\"n\">svd<\/span><span class=\"p\">(<\/span><span class=\"n\">M<\/span><span class=\"p\">)[<\/span><span class=\"mi\">1<\/span><span class=\"p\">]<\/span>\n    \n    <span class=\"n\">orig_M<\/span> <span class=\"o\">=<\/span> <span class=\"n\">M<\/span><span class=\"o\">.<\/span><span class=\"n\">copy<\/span><span class=\"p\">()<\/span>\n    <span class=\"n\">U<\/span><span class=\"p\">,<\/span> <span class=\"n\">S<\/span><span class=\"p\">,<\/span> <span class=\"n\">Vt<\/span> <span class=\"o\">=<\/span> <span class=\"n\">clear_svd<\/span><span class=\"p\">(<\/span><span class=\"n\">M<\/span><span class=\"p\">)<\/span>\n    \n    <span class=\"k\">return<\/span> <span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">allclose<\/span><span class=\"p\">(<\/span><span class=\"nb\">sorted<\/span><span class=\"p\">(<\/span><span class=\"n\">expected_S<\/span><span class=\"p\">),<\/span> <span class=\"nb\">sorted<\/span><span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">diag<\/span><span class=\"p\">(<\/span><span class=\"n\">S<\/span><span class=\"p\">)))<\/span> <span class=\"ow\">and<\/span>\n            <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">allclose<\/span><span class=\"p\">(<\/span><span class=\"n\">orig_M<\/span><span class=\"p\">,<\/span> <span class=\"n\">U<\/span><span class=\"o\">.<\/span><span class=\"n\">dot<\/span><span class=\"p\">(<\/span><span class=\"n\">S<\/span><span class=\"p\">)<\/span><span class=\"o\">.<\/span><span class=\"n\">dot<\/span><span class=\"p\">(<\/span><span class=\"n\">Vt<\/span><span class=\"p\">)))<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing code_cell rendered\">\n<div class=\"input\">\n<div class=\"prompt input_prompt\">In&nbsp;[13]:<\/div>\n<div class=\"inner_cell\">\n<div class=\"input_area\">\n<div class=\" highlight hl-ipython2\">\n<pre><span class=\"n\">test_matrices<\/span> <span class=\"o\">=<\/span> <span class=\"p\">[<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">diag<\/span><span class=\"p\">(<\/span><span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">arange<\/span><span class=\"p\">(<\/span><span class=\"mf\">5.0<\/span><span class=\"p\">)),<\/span>\n                 <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">arange<\/span><span class=\"p\">(<\/span><span class=\"mf\">16.0<\/span><span class=\"p\">)<\/span><span class=\"o\">.<\/span><span class=\"n\">reshape<\/span><span class=\"p\">(<\/span><span class=\"mi\">4<\/span><span class=\"p\">,<\/span><span class=\"mi\">4<\/span><span class=\"p\">),<\/span>\n\n                 <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">array<\/span><span class=\"p\">([[<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span> <span class=\"o\">-<\/span><span class=\"mf\">2.0<\/span><span class=\"p\">],<\/span>    <span class=\"c1\"># caused grief because of error in BBM\/Givens<\/span>\n                           <span class=\"p\">[<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span>    <span class=\"mi\">0<\/span><span class=\"p\">]]),<\/span>\n\n                 <span class=\"n\">np<\/span><span class=\"o\">.<\/span><span class=\"n\">array<\/span><span class=\"p\">([[<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">],<\/span>    <span class=\"c1\"># comes from applying 2 givens rotations to diag<\/span>\n                           <span class=\"p\">[<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mf\">2.12<\/span><span class=\"p\">,<\/span> <span class=\"o\">-<\/span><span class=\"mf\">2.12<\/span><span class=\"p\">],<\/span>\n                           <span class=\"p\">[<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mf\">1.414<\/span><span class=\"p\">,<\/span> <span class=\"mf\">1.414<\/span><span class=\"p\">],<\/span>\n                           <span class=\"p\">[<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">,<\/span><span class=\"mi\">0<\/span><span class=\"p\">]])]<\/span>\n\n<span class=\"k\">for<\/span> <span class=\"n\">A<\/span> <span class=\"ow\">in<\/span> <span class=\"n\">test_matrices<\/span><span class=\"p\">:<\/span>\n    <span class=\"k\">print<\/span> <span class=\"n\">test_check<\/span><span class=\"p\">(<\/span><span class=\"n\">A<\/span><span class=\"p\">)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"output_wrapper\">\n<div class=\"output\">\n<div class=\"output_area\">\n<div class=\"prompt\"><\/div>\n<div class=\"output_subarea output_stream output_stdout output_text\">\n<pre>True\nTrue\nTrue\nTrue\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\">\n<\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<h3>\nAdditional Resources<br \/>\n<\/h3>\n<p>\nYou can grab a <a href=\"http:\/\/drsfenner.org\/public\/notebooks\/ImplementingSVD.ipynb\">copy of this notebook<\/a>.\n<\/p>\n<p>\nEven better, you can <a href=\"http:\/\/nbviewer.ipython.org\/url\/drsfenner.org\/public\/notebooks\/ImplementingSVD.ipynb\">view it using nbviewer<\/a>.\n<\/p>\n<h3>\nLicense<br \/>\n<\/h3>\n<p>\nUnless otherwise noted, the contents of this notebook are under the following license. The code in the notebook should be considered part of the text (i.e., licensed and treated as as follows).\n<\/p>\n<p><a rel=\"license\" href=\"http:\/\/creativecommons.org\/licenses\/by-nc-sa\/4.0\/\"><img decoding=\"async\" alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https:\/\/i.creativecommons.org\/l\/by-nc-sa\/4.0\/88x31.png\" \/><\/a><br \/><span xmlns:dct=\"http:\/\/purl.org\/dc\/terms\/\" property=\"dct:title\">DrsFenner.org Blog And Notebooks<\/span> by <a xmlns:cc=\"http:\/\/creativecommons.org\/ns#\" href=\"drsfenner.org\" property=\"cc:attributionName\" rel=\"cc:attributionURL\">Mark and Barbara Fenner<\/a> is licensed under a <a rel=\"license\" href=\"http:\/\/creativecommons.org\/licenses\/by-nc-sa\/4.0\/\">Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License<\/a>.<br \/>Permissions beyond the scope of this license may be available at <a xmlns:cc=\"http:\/\/creativecommons.org\/ns#\" href=\"drsfenner.org\/blog\/about-and-contacts\" rel=\"cc:morePermissions\">drsfenner.org\/blog\/about-and-contacts<\/a>.<\/p>\n<\/div>\n<\/div>\n<\/div>\n","protected":false},"excerpt":{"rendered":"<p>At long last, I&#8217;ve gotten to my original goal: writing a relatively easy to understand version of your plain, vanilla SVD computation. I&#8217;m going to spend just a few minutes glossing (literally) over the theory behind the code and then dive into some implementation that builds on the previous posts. Thanks for reading!<\/p>\n","protected":false},"author":3,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[6,7],"tags":[],"class_list":["post-450","post","type-post","status-publish","format-standard","hentry","category-mrdr","category-sci-math-stat-python"],"_links":{"self":[{"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/450","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/comments?post=450"}],"version-history":[{"count":1,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/450\/revisions"}],"predecessor-version":[{"id":451,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/450\/revisions\/451"}],"wp:attachment":[{"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/media?parent=450"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/categories?post=450"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/tags?post=450"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}