{"id":436,"date":"2016-03-12T16:20:33","date_gmt":"2016-03-12T16:20:33","guid":{"rendered":"http:\/\/drsfenner.org\/blog\/?p=436"},"modified":"2016-03-12T16:27:15","modified_gmt":"2016-03-12T16:27:15","slug":"givens-rotations-and-the-case-of-the-blemished-bidiagonal-matrix","status":"publish","type":"post","link":"https:\/\/drsfenner.org\/blog\/2016\/03\/givens-rotations-and-the-case-of-the-blemished-bidiagonal-matrix\/","title":{"rendered":"Givens Rotations and the Case of the Blemished Bidiagonal Matrix"},"content":{"rendered":"
\n
\n<\/div>\n
\n
\n

Last time, we looked at using Givens rotations to perform a QR factorization of a matrix. Today, we’re going to be a little more selective and use Givens rotations to walk values off the edge of a special class of matrices <\/p>\n<\/div>\n<\/div>\n<\/div>\n

\n
\n<\/div>\n
\n
\n

Old Friends<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n

Here are few friends that we introduced last time<\/a>. I updated the interfaces a little bit: the parts that define the Givens rotation are tupled together and the target matrix is the left or right argument, depending on which multiplication we are doing.<\/p>\n<\/div>\n<\/div>\n<\/div>\n

\n
\n
In [1]:<\/div>\n
\n
\n
\n
import<\/span> numpy<\/span> as<\/span> np<\/span>\nimport<\/span> numpy.linalg<\/span> as<\/span> nla<\/span>\nimport<\/span> matplotlib.pyplot<\/span> as<\/span> plt<\/span>\n%<\/span>matplotlib<\/span> inline\nnp<\/span>.<\/span>set_printoptions<\/span>(<\/span>suppress<\/span>=<\/span>True<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [2]:<\/div>\n
\n
\n
\n
# GvL pg. 216 : algo 5.1.3<\/span>\ndef<\/span> zeroing_givens_coeffs<\/span>(<\/span>x<\/span>,<\/span>z<\/span>):<\/span>\n    ''' for the values x,z compute cos th, sin th <\/span>\n        s.t. applying a Givens rotation G(cos th,sin th) <\/span>\n             on 2 rows(or cols) with values x,z will<\/span>\n             maps x --> r and z --> 0'''<\/span>\n    if<\/span> z<\/span> ==<\/span> 0.0<\/span>:<\/span> # better to check for "small": abs(z) < np.finfo(np.double).eps<\/span>\n        return<\/span> 1.0<\/span>,<\/span>0.0<\/span>\n    r<\/span> =<\/span> np<\/span>.<\/span>hypot<\/span>(<\/span>x<\/span>,<\/span>z<\/span>)<\/span> # C99 hypot is safe for under\/overflow<\/span>\n    return<\/span> x<\/span>\/<\/span>r<\/span>,<\/span> -<\/span>z<\/span>\/<\/span>r<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [3]:<\/div>\n
\n
\n
\n
# GvL, pg. 216 .... Section 5.1.9<\/span>\ndef<\/span> left_givensT<\/span>((<\/span>c<\/span>,<\/span>s<\/span>,<\/span>r1<\/span>,<\/span>r2<\/span>),<\/span> A<\/span>):<\/span>\n    ''' update A <- G.T.dot(A) ... affects rows r1 and r2 '''<\/span>\n    givensT<\/span> =<\/span> np<\/span>.<\/span>array<\/span>([[<\/span> c<\/span>,<\/span> -<\/span>s<\/span>],<\/span>   # manually transposed <\/span>\n                        [<\/span> s<\/span>,<\/span>  c<\/span>]])<\/span>\n    A<\/span>[[<\/span>r1<\/span>,<\/span>r2<\/span>],:]<\/span> =<\/span> np<\/span>.<\/span>dot<\/span>(<\/span>givensT<\/span>,<\/span> A<\/span>[[<\/span>r1<\/span>,<\/span>r2<\/span>],:])<\/span>\n\ndef<\/span> right_givens<\/span>(<\/span>A<\/span>,<\/span> (<\/span>c<\/span>,<\/span>s<\/span>,<\/span>c1<\/span>,<\/span>c2<\/span>)):<\/span>\n    ''' update A <- A.dot(G) ... affects cols c1 and c2 '''<\/span>\n    givens<\/span> =<\/span> np<\/span>.<\/span>array<\/span>([[<\/span> c<\/span>,<\/span> s<\/span>],<\/span>\n                       [<\/span>-<\/span>s<\/span>,<\/span> c<\/span>]])<\/span>\n    A<\/span>[:,[<\/span>c1<\/span>,<\/span> c2<\/span>]]<\/span> =<\/span> np<\/span>.<\/span>dot<\/span>(<\/span>A<\/span>[:,[<\/span>c1<\/span>,<\/span> c2<\/span>]],<\/span> givens<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Zeroing in a Bidiagonal Matrix (Take 1)<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
\\[\\renewcommand{\\vec}[1]{\\mathbf{#1}}\\]<\/span> With the tools above and the lessons from the Givens QR post<\/a> we can come up with other uses for zeroing Givens rotations which I’ll abbreviate as \\(\\vec{G_0}(c,s,i,z)\\)<\/span>.<\/p>\n
We’re going to make use of \\(\\vec{G_0}\\)<\/span> in at least two different ways, when we perform an SVD (singular value decomposition). When we perform an SVD, we’ll take an intermediate step of making an (upper) bidiagonal matrix<\/em> — that is a matrix where the only non-zero entries are on the main and super-diagonals. Along the way, the clean bidiagonal matrix<\/em> may become blemished<\/em> and have a non-zero entry off the main upper bidiagonal band creating a blemished bidiagonal matrix<\/em> (BBM). We’ll use Givens rotations to "walk" the non-zero blemish out of the matrix. We can do that in three directions: across a row, up a column, or zig-zagging the main diagonal band.<\/p>\n
To get started, we need to easily create upper bidiagonal matrices.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [4]:<\/div>\n
\n
\n
\n
def<\/span> make_upper_bidiag<\/span>(<\/span>d1<\/span>,<\/span>d2<\/span>):<\/span>\n    ''' from two arrays, fill the diagonal and superdiagonal'''<\/span>\n    B<\/span> =<\/span> np<\/span>.<\/span>diag<\/span>(<\/span>d1<\/span>)<\/span>\n    B<\/span>.<\/span>flat<\/span>[<\/span>1<\/span>::<\/span>d1<\/span>.<\/span>size<\/span>+<\/span>1<\/span>]<\/span> =<\/span> d2<\/span> # awkward, but np.diag_indices is main diag only<\/span>\n    return<\/span> B<\/span>\n\ndef<\/span> extract_upper_bidiag<\/span>(<\/span>M<\/span>):<\/span>\n    '''from M, pull out the diagonal and superdiagonal <\/span>\n       like np.triu or np.tril would'''<\/span>\n    B<\/span> =<\/span> np<\/span>.<\/span>triu<\/span>(<\/span>M<\/span>)<\/span> *<\/span> np<\/span>.<\/span>tril<\/span>(<\/span>np<\/span>.<\/span>ones_like<\/span>(<\/span>M<\/span>),<\/span> 1<\/span>)<\/span>\n    return<\/span> B<\/span>\n    \nextract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>36.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>6<\/span>,<\/span>6<\/span>))<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
Out[4]:<\/div>\n
\n
array([[  0.,   1.,   0.,   0.,   0.,   0.],\n       [  0.,   7.,   8.,   0.,   0.,   0.],\n       [  0.,   0.,  14.,  15.,   0.,   0.],\n       [  0.,   0.,   0.,  21.,  22.,   0.],\n       [  0.,   0.,   0.,   0.,  28.,  29.],\n       [  0.,   0.,   0.,   0.,   0.,  35.]])<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
We might start by trying to "swap" the 8 with the 0 to its right. Basically, we’d try to take \\(x=0 \\rightarrow 8\\)<\/span> and \\(z=8 \\rightarrow 0\\)<\/span>; that is, calling zeroing_givens_coeffs(0, 8)<\/code>). Before we see the results, what rows or columns are going to be modified, if our \\(i,j\\)<\/span> are \\(1,2\\)<\/span> and we use a right Givens rotation?<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [5]:<\/div>\n
\n
\n
\n
def<\/span> move_nonzero_using_nextright<\/span>(<\/span>M<\/span>,<\/span> row<\/span>,<\/span> col<\/span>):<\/span>\n    cos_th<\/span>,<\/span>sin_th<\/span> =<\/span> zeroing_givens_coeffs<\/span>(<\/span>M<\/span>[<\/span>row<\/span>,<\/span>col<\/span>+<\/span>1<\/span>],<\/span> M<\/span>[<\/span>row<\/span>,<\/span>col<\/span>])<\/span>\n    right_givens<\/span>(<\/span>M<\/span>,<\/span> (<\/span>cos_th<\/span>,<\/span>sin_th<\/span>,<\/span> col<\/span>+<\/span>1<\/span>,<\/span> col<\/span>))<\/span>\n    \nB<\/span> =<\/span> extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>36.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>6<\/span>,<\/span>6<\/span>))<\/span>\nmove_nonzero_using_nextright<\/span>(<\/span>B<\/span>,<\/span> 1<\/span>,<\/span> 2<\/span>)<\/span>\nprint<\/span> B<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
[[  0.   1.   0.   0.   0.   0.]\n [  0.   7.   0.   8.   0.   0.]\n [  0.   0. -15.  14.   0.   0.]\n [  0.   0. -21.   0.  22.   0.]\n [  0.   0.   0.   0.  28.  29.]\n [  0.   0.   0.   0.   0.  35.]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
The Givens rotation makes use of the next column from the row,col<\/code> values we passed to move_nonzero_using_nextright<\/code>. So, the rotation affects columns 2 and 3 and introduces non-zeros at B[1,3], B[3,2]<\/code>. We modified more values than we really wanted. In particular, when we affect B[3,2]<\/code> we are modifying below<\/em> the diagonal and making more blemishes. Further, we don’t really care that the 8<\/code> stays an 8<\/code> when it is "swapped". So, instead of "swapping" a zero in, let’s find a value that is non-zero (and can stay any<\/em> non-zero value). We need look no further than the main diagonal.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Zeroing in a Bidiagonal Matrix (Take 2)<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
When I saw "swapping", I what I really want to do is introduce a zero off the diagonal and compute a new associated \\(r\\)<\/span> on the diagonal<\/em>. So, let’s find a new position to "swap" with. Any of the values on the diagonal are non-zero and<\/em>, as long as they stay non-zero, they don’t change the bidiagonal form of the matrix. Another aside: we actually care that the bidiagonal part is fully filled with non-zeros, or you might say, dense<\/em> – although that’s a slight abuse of the term.<\/p>\n
Anywho. Let’s try to target a "swap" between a spot \\(row,col\\)<\/span> and an entry on one of the two diagonal entries it shares in its row and column. When we do this, we will preserve much of the bidagonal matrix structure, since in a bidiagonal matrix, we already have entries on the diagonal. However, we will introduce or propogate a blemish. More on that in second. Here’s some code that works by using the diagonal elements. A mental quick-check: for right and left Givens rotations, do we use the shared row or column diagonal elements for our "swapping"?<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [6]:<\/div>\n
\n
\n
\n
def<\/span> givens_zero_wdiag_left<\/span>(<\/span>A<\/span>,<\/span> row<\/span>,<\/span> col<\/span>):<\/span>\n    ''' update two whole rows and at col: A[row,col] <- 0; A[col,col] <- r<\/span>\n        moves value up\/down to diagonal ... within the shared column'''<\/span>\n    c<\/span>,<\/span>s<\/span> =<\/span> zeroing_givens_coeffs<\/span>(<\/span>A<\/span>[<\/span>col<\/span>,<\/span>col<\/span>],<\/span> A<\/span>[<\/span>row<\/span>,<\/span>col<\/span>])<\/span>\n    G<\/span> =<\/span> (<\/span>c<\/span>,<\/span>s<\/span>,<\/span>col<\/span>,<\/span>row<\/span>)<\/span>  # swapped order to zero *second* arg<\/span>\n    left_givensT<\/span>(<\/span>G<\/span>,<\/span> A<\/span>)<\/span> # G tells affected *rows* - one col gets 0<\/span>\n    return<\/span> G<\/span>\n\ndef<\/span> givens_zero_wdiag_right<\/span>(<\/span>A<\/span>,<\/span> row<\/span>,<\/span> col<\/span>):<\/span>\n    ''' update two whole cols and at row: A[row,col] <- 0; A[row,row] <- r<\/span>\n        moves value left\/right to diagonal ... within the shared row'''<\/span>\n    c<\/span>,<\/span>s<\/span> =<\/span> zeroing_givens_coeffs<\/span>(<\/span>A<\/span>[<\/span>row<\/span>,<\/span>row<\/span>],<\/span> A<\/span>[<\/span>row<\/span>,<\/span>col<\/span>])<\/span>\n    G<\/span> =<\/span> (<\/span>c<\/span>,<\/span>s<\/span>,<\/span>row<\/span>,<\/span>col<\/span>)<\/span>\n    right_givens<\/span>(<\/span>A<\/span>,<\/span> G<\/span>)<\/span> # G tells affected *cols* - one row gets 0<\/span>\n    return<\/span> G<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
I mentioned, in the last post, that I didn’t like the coupling of two function calls (creating the coefficients and applying them). In the cell directly above, we couple the two calls inside a higher level function. This doesn’t really fix the argument ordering issue or the that fact that using a givens_zero_wdiag_left<\/code> with a right_givens<\/code> (application) would be legal code with a non-sensical result.<\/p>\n
If we really wanted to "fix" the software engineering issue, and stay with a functional (non-object oriented) solution, we could bundle the correct application (left or right) function with the result of the tuple arguments — which is basically what a class constructor would do. I did write a small set of classes to encapsulate this further, but for now I want to focus on the Givens parts, not the software engineering parts. The classes were even more useful because shortly we’ll have a special case where we can write a custom multiplication routine and save work on the left and right Givens applications. This gives quite a nice use of inheriting the basics and overriding the special case.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Back to the code we just wrote. On the left, the application looks like:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [7]:<\/div>\n
\n
\n
\n
B<\/span> =<\/span> extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>36.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>6<\/span>,<\/span>6<\/span>))<\/span>\ngivens_zero_wdiag_left<\/span>(<\/span>B<\/span>,<\/span> 1<\/span>,<\/span> 2<\/span>)<\/span>\nprint<\/span> B<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
[[  0.           1.           0.           0.           0.           0.        ]\n [  0.           6.07770199   0.          -7.44208408   0.           0.        ]\n [  0.           3.47297257  16.1245155   13.02364713   0.           0.        ]\n [  0.           0.           0.          21.          22.           0.        ]\n [  0.           0.           0.           0.          28.          29.        ]\n [  0.           0.           0.           0.           0.          35.        ]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
At this point, you’re ready to strangle me. That isn’t upper bidiagonal any more. But wait! What happens if we have and upper-bidiagonal with a zero on the diagonal?<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
In [8]:<\/div>\n
\n
\n
\n
B<\/span> =<\/span> extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>36.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>6<\/span>,<\/span>6<\/span>))<\/span>\nB<\/span>[<\/span>1<\/span>,<\/span>1<\/span>]<\/span> =<\/span> 0<\/span>\ngivens_zero_wdiag_left<\/span>(<\/span>B<\/span>,<\/span> 1<\/span>,<\/span> 2<\/span>)<\/span>\nprint<\/span> B<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n
[[  0.           1.           0.           0.           0.           0.        ]\n [  0.           0.           0.          -7.44208408   0.           0.        ]\n [  0.           0.          16.1245155   13.02364713   0.           0.        ]\n [  0.           0.           0.          21.          22.           0.        ]\n [  0.           0.           0.           0.          28.          29.        ]\n [  0.           0.           0.           0.           0.          35.        ]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
We move a value to the right and<\/em> keep the rest of the matrix upper bidiagonal (outside the blemish in the row we’re manipulating). That leads us to …<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Blemishes<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
I also left out one point about the form of the matrix we’ll be using attacking. In reality, le-subject-d’-attack<\/em> will be an upper-bidiagonal with a zero somewhere on the diagonal. The super-diagonal entry in the same row will be pushed across the row and out. Once it is off the super-diagonal, it is a "blemish" to the "pure" upper-bidiagonal matrix. Once we push it out, we’ll have two<\/em> upper-bidiagonal matrices (one above and one below) the row we modified.<\/p>\n
Here’s the general upper-bidiagonal matrix (no blemishes):<\/p>\n
\\[G=\\begin{bmatrix}
\nd_1    & s_1    & 0      &         & \\dots   & 0       \\\\
\n0      & d_2    & s_2    &         &         & 0       \\\\
\n       & \\ddots & \\ddots & \\ddots  & \\ddots  & \\vdots  \\\\
\n\\vdots &        &        & d_{n-2} & s_{n-2} & 0       \\\\
\n       &        &        &         & d_{n-1} & s_{n-1} \\\\
\n0      &        & \\dots  &         & 0       & d_n
\n\\end{bmatrix}\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
If we introduce a \\(0\\)<\/span> on the diagonal somewhere, we can talk about "walking" the superdiagonal element in that same row off to the right. Here, the two red values<\/font> are passed to zeroing_givens_coeffs<\/code> to compute the coefficient values that will zero the value in the upper row. To the right of an arrow, the blue elements<\/font> are newly computed updates. The \\(\\color{red}{B}\\)<\/span>, our blemish, has also been updated: \\(0 \\rightarrow \\color{red}{B}\\)<\/span>. Note, the \\(x\\)<\/span> and \\(B\\)<\/span> names merely indicate the zeroness of a value – two different \\(x\\)<\/span> or \\(B\\)<\/span> will have different values. Generally speaking \\(B\\)<\/span> will be different in different steps.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
\\[
\n\\renewcommand{\\red}[1]{\\color{red}{#1}}
\n\\renewcommand{\\blue}[1]{\\color{blue}{#1}}
\n\\begin{align}
\nA = \\begin{bmatrix}
\nx & x & 0 & 0 & 0 \\\\
\n0 & 0 & \\red{x} & 0 & 0 \\\\
\n0 & 0 & \\red{x} & x & 0 \\\\
\n0 & 0 & 0 & x & x \\\\
\n0 & 0 & 0 & 0 & x \\\\
\n\\end{bmatrix}
\n& \\xrightarrow{G^T_z(c,s,2,1) A} &
\n\\begin{bmatrix}
\nx & x & 0 & 0 & 0 \\\\
\n0 & 0 & \\blue{0} & \\red{B} & 0 \\\\
\n0 & 0 & \\blue{x} & \\blue{x} & 0 \\\\
\n0 & 0 & 0 & \\color{red}{x} & x \\\\
\n0 & 0 & 0 & 0 & x \\\\
\n\\end{bmatrix} \\\\
\n& \\xrightarrow{G^T_z(c,s,3,2) A} &
\n\\begin{bmatrix}
\nx & x & 0 & 0 & 0 \\\\
\n0 & 0 & 0 & \\blue{0} & \\red{B} \\\\
\n0 & 0 & x & x & 0 \\\\
\n0 & 0 & 0 & \\blue{x} & \\blue{x} \\\\
\n0 & 0 & 0 & 0 & \\color{red}{x} \\\\
\n\\end{bmatrix} \\\\
\n& \\xrightarrow{G^T_z(c,s,4,2) A} &
\n\\begin{bmatrix}
\nx & x & 0 & 0 & 0 \\\\
\n0 & 0 & 0 & 0 & \\blue{0} \\\\
\n0 & 0 & x & x & 0 \\\\
\n0 & 0 & 0 & x & x \\\\
\n0 & 0 & 0 & 0 & \\blue{x} \\\\
\n\\end{bmatrix}
\n\\end{align}
\n\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
When we move A[1,2]<\/code> to the right, we have to involve the \\((1,2)\\)<\/span>-plane. That is, we need a Givens rotation on the \\((1,2)\\)<\/span>-plane (and zeroing in the 1 plane). Since we’re applying a left-Givens, we are only updating two rows: 1<\/code> and 2<\/code>. Within those rows, the only nonzero entries are in columns 2<\/code> and 3<\/code>. So, using a similar argument to that we used for the QR algorithm last time, we can see that columns 0<\/code>, 1<\/code>, and 4<\/code> are unchanged and remain zero.<\/p>\n
We can get a bit more specific about the values for columns 2<\/code> and 3<\/code>. Let’s name the relevant \\(x\\)<\/span> values:<\/p>\n
    \n
  • \\(b\\)<\/span> the blemish. In \\(\\vec{A}\\)<\/span> it’s the only non-zero in row \\(1\\)<\/span>; it’s also the upper red \\(x\\)<\/span>.<\/li>\n
  • \\(d\\)<\/span> the diagonal value below<\/em> \\(b\\)<\/span>. In \\(\\vec{A}\\)<\/span>, the lower red \\(x\\)<\/span>.<\/li>\n
  • \\(f\\)<\/span> the friend. The value to the right of \\(d\\)<\/span> – it’s the super-diagonal super-Friend of \\(d\\)<\/span>. In \\(\\vec{A}\\)<\/span>, it is at A[2,3]<\/code>.<\/li>\n<\/ul>\n
    These are the only values from \\(\\vec{A}\\)<\/span> we need to calculate \\(\\vec{G^TA}\\)<\/span>. Everything else disappears to zero. In the last iteration, there is no value to the right of \\(d\\)<\/span>, so we can ignore it. We’ll can use this knowledge to optimize the matrix multiplication from two dot-products (the length of the whole rows) down to a couple of multiplications and additions.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Efficient Blemish Updating<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Based on the patterns of blue and red above, you might be thinking "gee, I don’t think we need the whole row" to do that update. We don’t. Turns out, we can do a "custom" multiplication routine for a Givens rotation against a blemished upper bidiagonal matrix<\/em>, or BBM. Here’s the same left Givens multiplication (as the first \\(\\vec{G^T_z}\\)<\/span> above) written out by hand but not doing any unnecessary computations with zeros. Remember, friend<\/code> is the resulting blue, non-zero superdiagonal entry involved in that step.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [9]:<\/div>\n
    \n
    \n
    \n
    B<\/span> =<\/span> extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>36.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>6<\/span>,<\/span>6<\/span>))<\/span>\nB<\/span>[<\/span>1<\/span>,<\/span>1<\/span>]<\/span> =<\/span> 0<\/span>\n\n# the setup:<\/span>\n# . . . . . <\/span>\n# . 0 B 0 0<\/span>\n# . 0 D F 0<\/span>\n# <\/span>\nrow<\/span>,<\/span>col<\/span> =<\/span> 1<\/span>,<\/span>2<\/span>\nblemish_loc<\/span> =<\/span> (<\/span>row<\/span>,<\/span>col<\/span>)<\/span>\ndiag_loc<\/span>    =<\/span> (<\/span>col<\/span>,<\/span>col<\/span>)<\/span>\nfriend_loc<\/span>  =<\/span> (<\/span>col<\/span>,<\/span> col<\/span>+<\/span>1<\/span>)<\/span>\n\nblemish<\/span> =<\/span> B<\/span>[<\/span>blemish_loc<\/span>]<\/span>\ndiag<\/span>    =<\/span> B<\/span>[<\/span>diag_loc<\/span>]<\/span>\nfriend<\/span>  =<\/span> B<\/span>[<\/span>friend_loc<\/span>]<\/span>\n\nc<\/span>,<\/span>s<\/span> =<\/span> zeroing_givens_coeffs<\/span>(<\/span>diag<\/span>,<\/span> blemish<\/span>)<\/span>\n\n# G^T = <\/span>\n#       [s  c]  --> zeroing row<\/span>\nB<\/span>[<\/span>blemish_loc<\/span>]<\/span> =<\/span> diag<\/span> *<\/span> s<\/span> +<\/span> blemish<\/span> *<\/span> c<\/span> # == 0<\/span>\nB<\/span>[<\/span>diag_loc<\/span>]<\/span>    =<\/span> diag<\/span> *<\/span> c<\/span> -<\/span> blemish<\/span> *<\/span> s<\/span>\n\nnew_blemish_loc<\/span> =<\/span> (<\/span>row<\/span>,<\/span> col<\/span>+<\/span>1<\/span>)<\/span>\nB<\/span>[<\/span>new_blemish_loc<\/span>]<\/span> =<\/span> s<\/span> *<\/span> friend<\/span>\nB<\/span>[<\/span>friend_loc<\/span>]<\/span>      =<\/span> c<\/span> *<\/span> friend<\/span>\n\nprint<\/span> B<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[  0.           1.           0.           0.           0.           0.        ]\n [  0.           0.           0.          -7.44208408   0.           0.        ]\n [  0.           0.          16.1245155   13.02364713   0.           0.        ]\n [  0.           0.           0.          21.          22.           0.        ]\n [  0.           0.           0.           0.          28.          29.        ]\n [  0.           0.           0.           0.           0.          35.        ]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Do we get the same thing using our machinery? We got the same result in In[8]<\/code>. We don’t need to be nearly as explicit as I was in In[9]<\/code>, nor do we need to save all of those values in temporaries. Trimming it down, we can do this:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [10]:<\/div>\n
    \n
    \n
    \n
    def<\/span> givens_zero_wdiag_left_on_BBM<\/span>(<\/span>BBM<\/span>,<\/span> row<\/span>,<\/span> col<\/span>):<\/span>\n    #assert really BBM<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> BBM<\/span>.<\/span>shape<\/span>\n    c<\/span>,<\/span>s<\/span> =<\/span> zeroing_givens_coeffs<\/span>(<\/span>BBM<\/span>[<\/span>col<\/span>,<\/span>col<\/span>],<\/span> BBM<\/span>[<\/span>row<\/span>,<\/span>col<\/span>])<\/span>\n\n    #       col   col+1<\/span>\n    # row   blem new_blem      G^T = <\/span>\n    #        <\/span>\n    # col   diag friend              [s  c] ---> zeroing row<\/span>\n\n    old_blemish<\/span> =<\/span> BBM<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span>\n    BBM<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> =<\/span> 0<\/span>                                  # from zeroing row<\/span>\n    BBM<\/span>[<\/span>col<\/span>,<\/span>col<\/span>]<\/span> =<\/span> BBM<\/span>[<\/span>col<\/span>,<\/span>col<\/span>]<\/span> *<\/span> c<\/span> -<\/span> old_blemish<\/span> *<\/span> s<\/span> # from top row<\/span>\n\n    if<\/span> col<\/span> <<\/span> n<\/span>-<\/span>1<\/span>:<\/span> # update these in all but last column<\/span>\n        BBM<\/span>[<\/span>row<\/span>,<\/span> col<\/span>+<\/span>1<\/span>]<\/span> =<\/span> s<\/span> *<\/span> BBM<\/span>[<\/span>col<\/span>,<\/span> col<\/span>+<\/span>1<\/span>]<\/span>\n        BBM<\/span>[<\/span>col<\/span>,<\/span> col<\/span>+<\/span>1<\/span>]<\/span> *=<\/span> c<\/span>\n    return<\/span> (<\/span>c<\/span>,<\/span>s<\/span>,<\/span>row<\/span>,<\/span>col<\/span>)<\/span> # == G<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [11]:<\/div>\n
    \n
    \n
    \n
    B<\/span> =<\/span> extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>sqrt<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>1.0<\/span>,<\/span>26.0<\/span>))<\/span>.<\/span>reshape<\/span>(<\/span>5<\/span>,<\/span>5<\/span>))<\/span>\nB<\/span>[<\/span>1<\/span>,<\/span>1<\/span>]<\/span> =<\/span> 0.0<\/span>\nprint<\/span> B<\/span>\n\ngivens_zero_wdiag_left_on_BBM<\/span>(<\/span>B<\/span>,<\/span> 1<\/span>,<\/span>2<\/span>)<\/span>\nprint<\/span> B<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[ 1.          1.41421356  0.          0.          0.        ]\n [ 0.          0.          2.82842712  0.          0.        ]\n [ 0.          0.          3.60555128  3.74165739  0.        ]\n [ 0.          0.          0.          4.35889894  4.47213595]\n [ 0.          0.          0.          0.          5.        ]]\n[[ 1.          1.41421356  0.          0.          0.        ]\n [ 0.          0.          0.         -2.30940108  0.        ]\n [ 0.          0.          4.58257569  2.94392029  0.        ]\n [ 0.          0.          0.          4.35889894  4.47213595]\n [ 0.          0.          0.          0.          5.        ]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    And finally, we can put that in a loop to walk a value all the way off the right edge of the matrix:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [12]:<\/div>\n
    \n
    \n
    \n
    def<\/span> walk_blemish_out_right<\/span>(<\/span>B<\/span>,<\/span> row<\/span>,<\/span> start_col<\/span>):<\/span> # usually start_col = row+1<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> B<\/span>.<\/span>shape<\/span>\n    for<\/span> col<\/span> in<\/span> xrange<\/span>(<\/span>start_col<\/span>,<\/span> n<\/span>):<\/span>\n        givens_zero_wdiag_left_on_BBM<\/span>(<\/span>B<\/span>,<\/span> row<\/span>,<\/span> col<\/span>)<\/span>\n\nB<\/span> =<\/span> extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>25.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>5<\/span>,<\/span>5<\/span>))<\/span>\nB<\/span>[<\/span>1<\/span>,<\/span>1<\/span>]<\/span> =<\/span> 0.0<\/span>\nprint<\/span> B<\/span>\n\nwalk_blemish_out_right<\/span>(<\/span>B<\/span>,<\/span> 1<\/span>,<\/span> 2<\/span>)<\/span>\nprint<\/span> B<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[  0.   1.   0.   0.   0.]\n [  0.   0.   7.   0.   0.]\n [  0.   0.  12.  13.   0.]\n [  0.   0.   0.  18.  19.]\n [  0.   0.   0.   0.  24.]]\n[[  0.           1.           0.           0.           0.        ]\n [  0.           0.           0.           0.           0.        ]\n [  0.           0.          13.89244399  11.22912571   0.        ]\n [  0.           0.           0.          19.15480973  17.85452348]\n [  0.           0.           0.           0.          24.86394963]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Woot! Err, well, at least it’s done. When you look at the result, you might notice that we took an "almost-full" bidiagonal matrix (with one bad row) and converted it into two "full" bidiagonal matrices separated by an empty row. If that looks like decomposing a larger problem into smaller ones, you’re right! Although, we are not leading up to the algorithm that is traditionally called "divide-and-conquer SVD".<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Walking Up<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    We have one other case to consider. If the zero on the diagonal is in the bottom-right corner, than we can’t<\/em> swap to the diagonal in that column because we’d end up taking a zero to a non-zero. So instead, we will walk the blemish up<\/em> the column, swapping to the diagonal in the same row. From one perspective, this is the opposite what we did when walking to the right.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Here’s a version using the full Givens multiplication. We can use it to check the output of a customized "blemish-aware" version in a minute.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [13]:<\/div>\n
    \n
    \n
    \n
    B<\/span> =<\/span> np<\/span>.<\/span>sqrt<\/span>(<\/span>extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>9.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>3<\/span>,<\/span>3<\/span>)))<\/span>\nB<\/span>[<\/span>2<\/span>,<\/span>2<\/span>]<\/span> =<\/span> 0<\/span>\nprint<\/span> B<\/span>\ngivens_zero_wdiag_right<\/span>(<\/span>B<\/span>,<\/span> 1<\/span>,<\/span> 2<\/span>)<\/span>\nprint<\/span> B<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[ 0.          1.          0.        ]\n [ 0.          2.          2.23606798]\n [ 0.          0.          0.        ]]\n[[ 0.          0.66666667 -0.74535599]\n [ 0.          3.          0.        ]\n [ 0.          0.          0.        ]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [14]:<\/div>\n
    \n
    \n
    \n
    def<\/span> givens_zero_wdiag_right_on_BBM<\/span>(<\/span>BBM<\/span>,<\/span> row<\/span>,<\/span> col<\/span>):<\/span>\n    #assert really BBM<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> BBM<\/span>.<\/span>shape<\/span>\n    c<\/span>,<\/span>s<\/span> =<\/span> zeroing_givens_coeffs<\/span>(<\/span>BBM<\/span>[<\/span>row<\/span>,<\/span>row<\/span>],<\/span> BBM<\/span>[<\/span>row<\/span>,<\/span>col<\/span>])<\/span>\n\n    #         row     col                 zero column<\/span>\n    # row-1   friend  new_blem    G = <\/span>\n    #        <\/span>\n    # row     diag    blem            [-s  c] <\/span>\n\n    old_blemish<\/span> =<\/span> BBM<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span>\n    BBM<\/span>[<\/span>row<\/span>,<\/span>col<\/span>]<\/span> =<\/span> 0<\/span>                                  # from zeroing col<\/span>\n    BBM<\/span>[<\/span>row<\/span>,<\/span>row<\/span>]<\/span> =<\/span> BBM<\/span>[<\/span>row<\/span>,<\/span>row<\/span>]<\/span> *<\/span> c<\/span> -<\/span> old_blemish<\/span> *<\/span> s<\/span> # from left col<\/span>\n\n    if<\/span> row<\/span> ><\/span> 0<\/span>:<\/span> # update these in all but last row<\/span>\n        BBM<\/span>[<\/span>row<\/span>-<\/span>1<\/span>,<\/span> col<\/span>]<\/span> =<\/span> s<\/span> *<\/span> BBM<\/span>[<\/span>row<\/span>-<\/span>1<\/span>,<\/span> row<\/span>]<\/span>\n        BBM<\/span>[<\/span>row<\/span>-<\/span>1<\/span>,<\/span> row<\/span>]<\/span> *=<\/span> c<\/span>\n    return<\/span> (<\/span>c<\/span>,<\/span>s<\/span>,<\/span>row<\/span>,<\/span>col<\/span>)<\/span> # == G<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [15]:<\/div>\n
    \n
    \n
    \n
    B<\/span> =<\/span> np<\/span>.<\/span>sqrt<\/span>(<\/span>extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>9.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>3<\/span>,<\/span>3<\/span>)))<\/span>\nB<\/span>[<\/span>2<\/span>,<\/span>2<\/span>]<\/span> =<\/span> 0.0<\/span>\nprint<\/span> B<\/span>\n\ngivens_zero_wdiag_right_on_BBM<\/span>(<\/span>B<\/span>,<\/span>1<\/span>,<\/span>2<\/span>)<\/span>\nprint<\/span> B<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[ 0.          1.          0.        ]\n [ 0.          2.          2.23606798]\n [ 0.          0.          0.        ]]\n[[ 0.          0.66666667 -0.74535599]\n [ 0.          3.          0.        ]\n [ 0.          0.          0.        ]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [16]:<\/div>\n
    \n
    \n
    \n
    def<\/span> walk_blemish_out_up<\/span>(<\/span>B<\/span>,<\/span> start_row<\/span>,<\/span> col<\/span>):<\/span> # usually start_row = col-1<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> B<\/span>.<\/span>shape<\/span>\n    for<\/span> row<\/span> in<\/span> reversed<\/span>(<\/span>xrange<\/span>(<\/span>start_row<\/span>+<\/span>1<\/span>)):<\/span>\n        givens_zero_wdiag_right_on_BBM<\/span>(<\/span>B<\/span>,<\/span> row<\/span>,<\/span> col<\/span>)<\/span>\n\nB<\/span> =<\/span> extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>25.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>5<\/span>,<\/span>5<\/span>))<\/span>\nB<\/span>[<\/span>4<\/span>,<\/span>4<\/span>]<\/span> =<\/span> 0.0<\/span>\nprint<\/span> B<\/span>\n\nwalk_blemish_out_up<\/span>(<\/span>B<\/span>,<\/span> 3<\/span>,<\/span> 4<\/span>)<\/span>\nprint<\/span> B<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[  0.   1.   0.   0.   0.]\n [  0.   6.   7.   0.   0.]\n [  0.   0.  12.  13.   0.]\n [  0.   0.   0.  18.  19.]\n [  0.   0.   0.   0.   0.]]\n[[  0.5849499    0.81106943   0.           0.           0.        ]\n [  0.           7.39764046   5.50226459   0.           0.        ]\n [  0.           0.          15.26644142   8.94068042   0.        ]\n [  0.           0.           0.          26.17250466   0.        ]\n [  0.           0.           0.           0.           0.        ]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Here, we’ve reduced the larger matrix to a slightly smaller matrix. This allows us to solve a single, slightly smaller, problem. It’s just one subproblem, instead of the two subproblems we got when we walked a blemish out along a row.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    And Finally, Cleaning up the Superdiagonal<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    In G&VL, they call the decomposing of a large problem into two smaller problems "decoupling". And, I’ll call the reduction part reduce. So, depending on where we find a zero on the main diagonal, we’ll either decouple or reduce. We can put those together:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [17]:<\/div>\n
    \n
    \n
    \n
    def<\/span> decouple_or_reduce<\/span>(<\/span>B<\/span>,<\/span> zero_idxs<\/span>):<\/span>\n    # see GvL, pg. 454 and pray for mojo<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> B<\/span>.<\/span>shape<\/span>\n    for<\/span> kk<\/span> in<\/span> zero_idxs<\/span>:<\/span>\n        if<\/span> kk<\/span> <<\/span> n<\/span>-<\/span>1<\/span>:<\/span>\n            walk_blemish_out_right<\/span>(<\/span>B<\/span>,<\/span> kk<\/span>,<\/span> kk<\/span>+<\/span>1<\/span>)<\/span> # decouple<\/span>\n        elif<\/span> kk<\/span> ==<\/span> n<\/span>-<\/span>1<\/span>:<\/span>\n            walk_blemish_out_up<\/span>(<\/span>B<\/span>,<\/span> kk<\/span>-<\/span>1<\/span>,<\/span> kk<\/span>)<\/span>    # reduce<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [18]:<\/div>\n
    \n
    \n
    \n
    B<\/span> =<\/span> extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>25.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>5<\/span>,<\/span>5<\/span>))<\/span>\nB<\/span>[<\/span>4<\/span>,<\/span>4<\/span>]<\/span> =<\/span> 0.0<\/span>\nprint<\/span> B<\/span>\ndecouple_or_reduce<\/span>(<\/span>B<\/span>,<\/span> [<\/span>4<\/span>])<\/span>\nprint<\/span> B<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[  0.   1.   0.   0.   0.]\n [  0.   6.   7.   0.   0.]\n [  0.   0.  12.  13.   0.]\n [  0.   0.   0.  18.  19.]\n [  0.   0.   0.   0.   0.]]\n[[  0.5849499    0.81106943   0.           0.           0.        ]\n [  0.           7.39764046   5.50226459   0.           0.        ]\n [  0.           0.          15.26644142   8.94068042   0.        ]\n [  0.           0.           0.          26.17250466   0.        ]\n [  0.           0.           0.           0.           0.        ]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [19]:<\/div>\n
    \n
    \n
    \n
    B<\/span> =<\/span> extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>25.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>5<\/span>,<\/span>5<\/span>))<\/span>\nB<\/span>[<\/span>1<\/span>,<\/span>1<\/span>]<\/span> =<\/span> 0.0<\/span>\nprint<\/span> B<\/span>\ndecouple_or_reduce<\/span>(<\/span>B<\/span>,<\/span> [<\/span>1<\/span>])<\/span>\nprint<\/span> B<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[  0.   1.   0.   0.   0.]\n [  0.   0.   7.   0.   0.]\n [  0.   0.  12.  13.   0.]\n [  0.   0.   0.  18.  19.]\n [  0.   0.   0.   0.  24.]]\n[[  0.           1.           0.           0.           0.        ]\n [  0.           0.           0.           0.           0.        ]\n [  0.           0.          13.89244399  11.22912571   0.        ]\n [  0.           0.           0.          19.15480973  17.85452348]\n [  0.           0.           0.           0.          24.86394963]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    One Other Walk Out<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    One other scenario we need to deal with (in computing the SVD), is bouncing or zig-zagging a blemish down the band of a "dense" bidiagonal matrix. The red values are used to compute the zeroing coefficients for that \\(G_z\\)<\/span>. Here’s what it looks like:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    \\[
    \n\\begin{align}
    \nA =
    \n\\begin{bmatrix}
    \n\\red{x} & x & 0 & 0 & 0 \\\\
    \n\\red{B} & x & x & 0 & 0 \\\\
    \n0 & 0 & x & x & 0 \\\\
    \n0 & 0 & 0 & x & x \\\\
    \n0 & 0 & 0 & 0 & x \\\\
    \n\\end{bmatrix}
    \n& \\xrightarrow{G^T_z(c,s,0,1) A} &
    \n\\begin{bmatrix}
    \nx & \\red{x} & \\red{B} & 0 & 0 \\\\
    \n0 & x & x & 0 & 0 \\\\
    \n0 & 0 & x & x & 0 \\\\
    \n0 & 0 & 0 & x & x \\\\
    \n0 & 0 & 0 & 0 & x \\\\
    \n\\end{bmatrix}
    \n& \\xrightarrow{A G_z(c,s,1,2)} &
    \n\\begin{bmatrix}
    \nx & x & 0 & 0 & 0 \\\\
    \n0 & \\red{x} & x & 0 & 0 \\\\
    \n0 & \\red{B} & x & x & 0 \\\\
    \n0 & 0 & 0 & x & x \\\\
    \n0 & 0 & 0 & 0 & x \\\\
    \n\\end{bmatrix} \\\\
    \n& \\xrightarrow{G^T_z(c,s,1,2) A} &
    \n\\begin{bmatrix}
    \nx & x & 0 & 0 & 0 \\\\
    \n0 & x & \\red{x} & \\red{B} & 0 \\\\
    \n0 & 0 & x & x & 0 \\\\
    \n0 & 0 & 0 & x & x \\\\
    \n0 & 0 & 0 & 0 & x \\\\
    \n\\end{bmatrix}
    \n& \\xrightarrow{A G_z(c,s,2,3)} &
    \n\\begin{bmatrix}
    \nx & x & 0 & 0 & 0 \\\\
    \n0 & x & x & 0 & 0 \\\\
    \n0 & x & \\red{x} & x & 0 \\\\
    \n0 & 0 & \\red{B} & x & x \\\\
    \n0 & 0 & 0 & 0 & x \\\\
    \n\\end{bmatrix} \\\\
    \n& \\xrightarrow{G^T_z(c,s,2,3) A} &
    \n\\begin{bmatrix}
    \nx & x & 0 & 0 & 0 \\\\
    \n0 & x & x & 0 & 0 \\\\
    \n0 & 0 & x & \\red{x} & \\red{B} \\\\
    \n0 & 0 & 0 & x & x \\\\
    \n0 & 0 & 0 & 0 & x \\\\
    \n\\end{bmatrix}
    \n& \\xrightarrow{A G_z(c,s,3,4)} &
    \n\\begin{bmatrix}
    \nx & x & 0 & 0 & 0 \\\\
    \n0 & x & x & 0 & 0 \\\\
    \n0 & x & x & x & 0 \\\\
    \n0 & 0 & 0 & \\red{x} & x \\\\
    \n0 & 0 & 0 & \\red{B} & x \\\\
    \n\\end{bmatrix}
    \n\\end{align}
    \n\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    There is a very important difference between this zig-zagging – it looks like we’ve got a drunken blemish on our hands – and the "straighline" walkouts we performed above. Instead of "swapping" to the diagonal, we are going to "swap" to an adjacent element on the band. The adjacent element will either be (1) in the same row and to the left of us or (2) in the same column and above us. If we don’t use those neighbors or friends, we end up with more blemishes. Since we have a BBM we use those friends, we keep one, and only one, blemish.<\/p>\n
    However, using those friend positions, means that we need some slightly different left and right Givens rotations. Here’s the corresponding functions that perform these two steps:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [20]:<\/div>\n
    \n
    \n
    \n
    def<\/span> givens_zero_adjacent_onband_left<\/span>(<\/span>A<\/span>,<\/span> row<\/span>,<\/span> col<\/span>):<\/span>\n    # affecting two rows: need the row above me<\/span>\n    # assert col==row-1 ... below diag ... row>col<\/span>\n    c<\/span>,<\/span>s<\/span> =<\/span> zeroing_givens_coeffs<\/span>(<\/span>A<\/span>[<\/span>row<\/span>-<\/span>1<\/span>,<\/span>col<\/span>],<\/span> A<\/span>[<\/span>row<\/span>,<\/span>col<\/span>])<\/span> # same col<\/span>\n    G<\/span> =<\/span> (<\/span>c<\/span>,<\/span>s<\/span>,<\/span>row<\/span>-<\/span>1<\/span>,<\/span>row<\/span>)<\/span>\n    left_givensT<\/span>(<\/span>G<\/span>,<\/span> A<\/span>)<\/span> # G tells affected *rows* - one col gets 0<\/span>\n    return<\/span> G<\/span>\n\ndef<\/span> givens_zero_adjacent_onband_right<\/span>(<\/span>A<\/span>,<\/span> row<\/span>,<\/span> col<\/span>):<\/span>\n    # affecting two cols:  need col left of me<\/span>\n    # assert col-2==row ... above diag ... col>row<\/span>\n    c<\/span>,<\/span>s<\/span> =<\/span> zeroing_givens_coeffs<\/span>(<\/span>A<\/span>[<\/span>row<\/span>,<\/span>col<\/span>-<\/span>1<\/span>],<\/span> A<\/span>[<\/span>row<\/span>,<\/span>col<\/span>])<\/span> # same row<\/span>\n    G<\/span> =<\/span> (<\/span>c<\/span>,<\/span>s<\/span>,<\/span>col<\/span>-<\/span>1<\/span>,<\/span>col<\/span>)<\/span>\n    right_givens<\/span>(<\/span>A<\/span>,<\/span> G<\/span>)<\/span> # G tells affected *cols* - one row gets 0<\/span>\n    return<\/span> G<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    As with givens_zero_wdiag_left<\/code> and givens_zero_wdiag_right<\/code>, when applied to a BBM these have a very efficient form. I’m going to leave that as an exercise to the highly motivated reader. But, it will be very close to the forms for givens_zero_wdiag_left_on_BBM<\/code> (and for the right side application). If you follow the Givens rotations under "One Other Walkout", that seems pretty clear. We can turn that into code pretty quickly (adding one last left rotation to clear out the blemish in the bottom row):<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [21]:<\/div>\n
    \n
    \n
    \n
    def<\/span> bounce_blemish_down_diagonal<\/span>(<\/span>BBM<\/span>):<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> BBM<\/span>.<\/span>shape<\/span>\n    for<\/span> k<\/span> in<\/span> xrange<\/span>(<\/span>n<\/span>-<\/span>2<\/span>):<\/span>\n        givens_zero_adjacent_onband_left<\/span>(<\/span> BBM<\/span>,<\/span> k<\/span>+<\/span>1<\/span>,<\/span> k<\/span>)<\/span>\n        givens_zero_adjacent_onband_right<\/span>(<\/span>BBM<\/span>,<\/span> k<\/span>,<\/span>   k<\/span>+<\/span>2<\/span>)<\/span>\n\n    givens_zero_adjacent_onband_left<\/span>(<\/span>BBM<\/span>,<\/span> n<\/span>-<\/span>1<\/span>,<\/span> n<\/span>-<\/span>2<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [22]:<\/div>\n
    \n
    \n
    \n
    A<\/span> =<\/span> extract_upper_bidiag<\/span>(<\/span>np<\/span>.<\/span>arange<\/span>(<\/span>1.0<\/span>,<\/span> 17.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>4<\/span>,<\/span>4<\/span>))<\/span>\nA<\/span>[<\/span>1<\/span>,<\/span>0<\/span>]<\/span> =<\/span> 99<\/span>\nprint<\/span> A<\/span>\nbounce_blemish_down_diagonal<\/span>(<\/span>A<\/span>)<\/span>\nprint<\/span> A<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[  1.   2.   0.   0.]\n [ 99.   6.   7.   0.]\n [  0.   0.  11.  12.]\n [  0.   0.   0.  16.]]\n[[ 99.00505038   9.23223353  -0.           0.        ]\n [  0.           8.42736916  13.72460309   0.        ]\n [  0.           0.          14.11692105   7.60672004]\n [  0.           0.           0.          -2.86896033]]\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    \nAdditional Resources
    \n<\/h3>\n
    \nYou can grab a copy of this notebook<\/a>.\n<\/p>\n
    \nEven better, you can view it using nbviewer<\/a>.\n<\/p>\n
    \nLicense
    \n<\/h3>\n
    \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
    \"Creative<\/a>
    DrsFenner.org Blog And Notebooks<\/span> by     Mark and Barbara Fenner<\/a> is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License<\/a>.
    Permissions beyond the scope of this license may be available at     drsfenner.org\/blog\/about-and-contacts<\/a>.<\/p>\n<\/div>\n<\/div>\n<\/div>\n","protected":false},"excerpt":{"rendered":"
    Last time, we looked at using Givens rotations to perform a QR factorization of a matrix. Today, we’re going to be a little more selective and use Givens rotations to walk values off the edge of a special class of matrices<\/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-436","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\/436","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=436"}],"version-history":[{"count":1,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/436\/revisions"}],"predecessor-version":[{"id":437,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/436\/revisions\/437"}],"wp:attachment":[{"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/media?parent=436"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/categories?post=436"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/tags?post=436"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}