{"id":445,"date":"2016-03-15T14:02:54","date_gmt":"2016-03-15T14:02:54","guid":{"rendered":"http:\/\/drsfenner.org\/blog\/?p=445"},"modified":"2016-03-15T14:03:59","modified_gmt":"2016-03-15T14:03:59","slug":"householder-bidiagonalization","status":"publish","type":"post","link":"http:\/\/drsfenner.org\/blog\/2016\/03\/householder-bidiagonalization\/","title":{"rendered":"Householder Bidiagonalization"},"content":{"rendered":"
\n
\n<\/div>\n
\n
\n

We’re going to make use of Householder reflections<\/a> to turn a generic matrix into a bidiagonal matrix<\/a>. I’ve laid a lot of foundation for these topics in the posts I just linked. Check them out! Onward.<\/p>\n<\/div>\n<\/div>\n<\/div>\n

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

Reducing a Generic Matrix to a Bidiagonal Matrix Using Householder Reflections<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n

Why would we want to do this? Well, it’s a long chain of logic (and this is abbreviated) but with an input matrix \\(A\\)<\/span> that we eventually want to compute the SVD of \\(A=U_S D V^T_S\\)<\/span>. BTW, I’m using \\(D\\)<\/span> instead of the usual \\(\\Sigma\\)<\/span> because I don’t want folks thinking it is a variance\/covariance matrix and it is<\/em> diagonal, in any case.<\/p>\n

    \n
  1. If we have \\(U^T_B A V_B = B\\)<\/span> with \\(B\\)<\/span> upper-bidiagonal and \\(U_B,V_B\\)<\/span> orthogonal, then \\(V^T_B(A^T A)V_B=B^T B=T\\)<\/span> with \\(T\\)<\/span> tridiagonal and symmetric. Amazingly (serendipity?), the \\(U^T\\)<\/span> disappers since: \\((XYZ)^T=X^T Y^T Z^T\\)<\/span> and \\(U\\)<\/span> is orthogonal (\\(U U^T = I\\)<\/span>).<\/li>\n
  2. The singular values of \\(A\\)<\/span> are the positive square roots of the non-zero eigenvalues of \\(A^TA\\)<\/span>. Take a look at the inner term of the second equation in point 1.<\/li>\n
  3. \\(A^TA\\)<\/span> is similar — which means they have the same eigenvalues — to \\(B^TB\\)<\/span> and \\(T\\)<\/span>.<\/li>\n
  4. A symmetric tridiagonal matrix gives rise to a particularly nice characteristic equation which allows "easy" determination of eigenvalues (they end up forming a Sturm sequence<\/em> which has nice properties).<\/li>\n<\/ol>\n

    I don’t want to leave a misleading impression. We won’t end up explicitly using \\(B^TB\\)<\/span>. The reasons are the same reasons we avoid that form in solving least squares using the direct normal equations: it has bad numerical properties and can be large. However, we do need the bidiagonal form \\(B\\)<\/span>. Here’s make_house_vec<\/code> from the post on using Householder reflections to perform QR<\/a><\/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>\nnp<\/span>.<\/span>set_printoptions<\/span>(<\/span>suppress<\/span>=<\/span>True<\/span>)<\/span>\n\n# G&vL 3rd pg. 210<\/span>\ndef<\/span> make_house_vec<\/span>(<\/span>x<\/span>):<\/span>\n    n<\/span> =<\/span> x<\/span>.<\/span>shape<\/span>[<\/span>0<\/span>]<\/span>\n    dot_1on<\/span> =<\/span> x<\/span>[<\/span>1<\/span>:]<\/span>.<\/span>dot<\/span>(<\/span>x<\/span>[<\/span>1<\/span>:])<\/span>\n\n    # v is our return vector; we hack on v[0]<\/span>\n    v<\/span> =<\/span> np<\/span>.<\/span>copy<\/span>(<\/span>x<\/span>)<\/span>\n    v<\/span>[<\/span>0<\/span>]<\/span> =<\/span> 1.0<\/span>\n    \n    if<\/span> dot_1on<\/span> <<\/span> np<\/span>.<\/span>finfo<\/span>(<\/span>float<\/span>)<\/span>.<\/span>eps<\/span>:<\/span>\n        beta<\/span> =<\/span> 0.0<\/span>\n    else<\/span>:<\/span>\n        # apply Parlett's formula (G&vL page 210) for safe v_0 = x_0 - norm(X) <\/span>\n        norm_x<\/span>=<\/span> np<\/span>.<\/span>sqrt<\/span>(<\/span>x<\/span>[<\/span>0<\/span>]<\/span>**<\/span>2<\/span> +<\/span> dot_1on<\/span>)<\/span>\n        if<\/span> x<\/span>[<\/span>0<\/span>]<\/span> <=<\/span> 0<\/span>:<\/span>\n            v<\/span>[<\/span>0<\/span>]<\/span> =<\/span> x<\/span>[<\/span>0<\/span>]<\/span> -<\/span> norm_x<\/span>\n        else<\/span>:<\/span>\n            v<\/span>[<\/span>0<\/span>]<\/span> =<\/span> -<\/span>dot_1on<\/span> \/<\/span> (<\/span>x<\/span>[<\/span>0<\/span>]<\/span> +<\/span> norm_x<\/span>)<\/span>\n        beta<\/span> =<\/span> 2<\/span> *<\/span> v<\/span>[<\/span>0<\/span>]<\/span>**<\/span>2<\/span> \/<\/span> (<\/span>dot_1on<\/span> +<\/span> v<\/span>[<\/span>0<\/span>]<\/span>**<\/span>2<\/span>)<\/span>\n        v<\/span> =<\/span> v<\/span> \/<\/span> v<\/span>[<\/span>0<\/span>]<\/span>\n    return<\/span> v<\/span>,<\/span> beta<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Longwinded and Explicit<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Using Householder vectors, we’re going to introduce "lots" of zeros into colums of our target matrix. This is very similar to the Householder QR process<\/a>. However, instead of simply zeroing out below<\/em> the diagonal one column at a time, we’re also<\/em> going to zero out above the superdiagonal. We’re also going to work down the column and across the row. So, for position \\((p,p)\\)<\/span> on the diagonal of a \\((m,n)\\)<\/span> matrix, we’ll zero out elements in that column (using zero based indexing):<\/p>\n
    \\[(p+1,p),\\ (p+2,p),\\ \\ldots\\ ,\\ (m-1,p)\\]<\/span><\/p>\n
    and<\/em> we’ll zero out elements in that row:<\/p>\n
    \\[(p,p+2),\\ (p,p+3),\\ \\ldots\\ ,\\ (p, n-1)\\]<\/span><\/p>\n
    We start across the row at \\(p+2\\)<\/span> because we need to skip the superdiagonal element. To keep the implementation a bit more direct, we’re going to hold off on one optimization. Instead of "packing" the reflection matrices back into A<\/code>, we’re going to compute them explicitly. You don’t need to do this in practice, but hold that thought.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [2]:<\/div>\n
    \n
    \n
    \n
    def<\/span> full_house<\/span>(<\/span>n<\/span>,<\/span> col<\/span>,<\/span> v<\/span>,<\/span> beta<\/span>):<\/span>\n    ''' for size n, apply a Householder vector v in the lower right corner of <\/span>\n        I_n to get a full-sized matrix with a smaller Householder matrix component'''<\/span>\n    full<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>n<\/span>)<\/span>\n    full<\/span>[<\/span>col<\/span>:,<\/span> col<\/span>:]<\/span> -=<\/span> beta<\/span> *<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>v<\/span>,<\/span>v<\/span>)<\/span>\n    return<\/span> full<\/span>\n\n# G&VL Algo. 5.4.2 with explicit reflections<\/span>\ndef<\/span> house_bidiag_explicit_UV<\/span>(<\/span>A<\/span>):<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>\n    assert<\/span> m<\/span> >=<\/span> n<\/span>\n    U<\/span>,<\/span>Vt<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>m<\/span>),<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>n<\/span>)<\/span>\n    \n    for<\/span> col<\/span> in<\/span> xrange<\/span>(<\/span>n<\/span>):<\/span>\n        v<\/span>,<\/span> beta<\/span> =<\/span> make_house_vec<\/span>(<\/span>A<\/span>[<\/span>col<\/span>:,<\/span>col<\/span>])<\/span>\n        A<\/span>[<\/span>col<\/span>:,<\/span>col<\/span>:]<\/span> =<\/span> (<\/span>np<\/span>.<\/span>eye<\/span>(<\/span>m<\/span>-<\/span>col<\/span>)<\/span> -<\/span> beta<\/span> *<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>v<\/span>,<\/span>v<\/span>))<\/span>.<\/span>dot<\/span>(<\/span>A<\/span>[<\/span>col<\/span>:,<\/span>col<\/span>:])<\/span>\n        Q<\/span> =<\/span> full_house<\/span>(<\/span>m<\/span>,<\/span> col<\/span>,<\/span> v<\/span>,<\/span> beta<\/span>)<\/span>\n        U<\/span> =<\/span> U<\/span>.<\/span>dot<\/span>(<\/span>Q<\/span>)<\/span>\n        \n        if<\/span> col<\/span> <=<\/span> n<\/span>-<\/span>2<\/span>:<\/span>\n            # transpose here, reflection for zeros above diagonal in A<\/span>\n            # col+1 keeps us off the super diagonal<\/span>\n            v<\/span>,<\/span>beta<\/span> =<\/span> make_house_vec<\/span>(<\/span>A<\/span>[<\/span>col<\/span>,<\/span>col<\/span>+<\/span>1<\/span>:]<\/span>.<\/span>T<\/span>)<\/span>\n            A<\/span>[<\/span>col<\/span>:,<\/span>col<\/span>+<\/span>1<\/span>:]<\/span> =<\/span> A<\/span>[<\/span>col<\/span>:,<\/span> col<\/span>+<\/span>1<\/span>:]<\/span>.<\/span>dot<\/span>(<\/span>np<\/span>.<\/span>eye<\/span>(<\/span>n<\/span>-<\/span>(<\/span>col<\/span>+<\/span>1<\/span>))<\/span> -<\/span> beta<\/span> *<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>v<\/span>,<\/span>v<\/span>))<\/span>\n            P<\/span> =<\/span> full_house<\/span>(<\/span>n<\/span>,<\/span> col<\/span>+<\/span>1<\/span>,<\/span> v<\/span>,<\/span> beta<\/span>)<\/span>\n            Vt<\/span> =<\/span> P<\/span>.<\/span>dot<\/span>(<\/span>Vt<\/span>)<\/span>\n    return<\/span> U<\/span>,<\/span> A<\/span>,<\/span> Vt<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [3]:<\/div>\n
    \n
    \n
    \n
    A<\/span> =<\/span> np<\/span>.<\/span>arange<\/span>(<\/span>1<\/span>,<\/span>13.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>4<\/span>,<\/span>3<\/span>)<\/span>\nA_test<\/span> =<\/span> np<\/span>.<\/span>copy<\/span>(<\/span>A<\/span>)<\/span>\nprint<\/span> A<\/span>\n\nU<\/span>,<\/span> B<\/span>,<\/span> Vt<\/span> =<\/span> house_bidiag_explicit_UV<\/span>(<\/span>A<\/span>)<\/span>\nprint<\/span> "B:<\/span>\\n<\/span>"<\/span>,<\/span> B<\/span>\nprint<\/span> "U:<\/span>\\n<\/span>"<\/span>,<\/span> U<\/span>\nprint<\/span> "Vt:<\/span>\\n<\/span>"<\/span>,<\/span> Vt<\/span>\n\nprint<\/span> "should equal input:"<\/span>,<\/span> np<\/span>.<\/span>allclose<\/span>(<\/span>U<\/span>.<\/span>dot<\/span>(<\/span>B<\/span>)<\/span>.<\/span>dot<\/span>(<\/span>Vt<\/span>),<\/span> A_test<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[  1.   2.   3.]\n [  4.   5.   6.]\n [  7.   8.   9.]\n [ 10.  11.  12.]]\nB:\n[[ 12.88409873  21.87643283   0.        ]\n [  0.           2.24623524  -0.61328133]\n [ -0.           0.           0.        ]\n [ -0.           0.           0.        ]]\nU:\n[[ 0.07761505  0.83305216 -0.54737648  0.01946767]\n [ 0.31046021  0.45123659  0.74434565  0.38203344]\n [ 0.54330537  0.06942101  0.15343813 -0.8224699 ]\n [ 0.77615053 -0.31239456 -0.35040731  0.42096879]]\nVt:\n[[ 1.          0.          0.        ]\n [ 0.          0.66700225  0.7450557 ]\n [ 0.          0.7450557  -0.66700225]]\nshould equal input: True\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Refactored and Explicit<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Reading (let alone coding) the row<\/code>, col<\/code>, +1<\/code>, n-col<\/code> values in that version is painful. It might be worthwhile, for performance purposes, to have so many long-winded indexing operations. Those operations are fast (-er than python level function calls), but lack human-readable semantics — I apologize to our numerical wizard siblings who speak in binary. So, here’s another look at it, with some functions to apply the Householder vectors to both the input matrix and the rotations that we are carrying along for the ride.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [4]:<\/div>\n
    \n
    \n
    \n
    def<\/span> apply_house_left<\/span>(<\/span>submatrix<\/span>,<\/span> hvec<\/span>,<\/span> beta<\/span>,<\/span> rotation<\/span>,<\/span> fullrows<\/span>):<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> submatrix<\/span>.<\/span>shape<\/span>\n    \n    # A needs a reduced H (just the bottom-right); <\/span>\n    # rotation needs a full H<\/span>\n    fullH<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>fullrows<\/span>)<\/span>\n    partH<\/span> =<\/span> fullH<\/span>[<\/span>fullrows<\/span>-<\/span>m<\/span>:,<\/span> fullrows<\/span>-<\/span>m<\/span>:]<\/span>\n    partH<\/span> -=<\/span>  beta<\/span>*<\/span>np<\/span>.<\/span>outer<\/span>(<\/span>hvec<\/span>,<\/span>hvec<\/span>)<\/span>\n    \n    # mutate instead of (new) object assignment<\/span>\n    submatrix<\/span>[:]<\/span> =<\/span> partH<\/span>.<\/span>dot<\/span>(<\/span>submatrix<\/span>)<\/span>\n    rotation<\/span>[:]<\/span> =<\/span> rotation<\/span>.<\/span>dot<\/span>(<\/span>fullH<\/span>)<\/span>\n\n\ndef<\/span> apply_house_right<\/span>(<\/span>submatrix<\/span>,<\/span> hvec<\/span>,<\/span> beta<\/span>,<\/span> rotation<\/span>,<\/span> fullcols<\/span>):<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> submatrix<\/span>.<\/span>shape<\/span>\n\n    # A needs a reduced H (just the bottom-right); <\/span>\n    # rotation needs a full H<\/span>\n    fullH<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>fullcols<\/span>)<\/span>\n    partH<\/span> =<\/span> fullH<\/span>[<\/span>fullcols<\/span>-<\/span>n<\/span>:,<\/span> fullcols<\/span>-<\/span>n<\/span>:]<\/span> # fillable = fullcols - m  <\/span>\n    partH<\/span> -=<\/span> beta<\/span>*<\/span>np<\/span>.<\/span>outer<\/span>(<\/span>hvec<\/span>,<\/span>hvec<\/span>)<\/span>\n    \n    # mutate instead of (new) object assignment<\/span>\n    submatrix<\/span>[:]<\/span> =<\/span> submatrix<\/span>.<\/span>dot<\/span>(<\/span>partH<\/span>)<\/span>\n    rotation<\/span>[:]<\/span>  =<\/span> fullH<\/span>.<\/span>dot<\/span>(<\/span>rotation<\/span>)<\/span>\n\ndef<\/span> hbd_simple<\/span>(<\/span>A<\/span>):<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>\n    U<\/span>,<\/span>Vt<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>m<\/span>),<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>n<\/span>)<\/span>\n    \n    for<\/span> col<\/span> in<\/span> xrange<\/span>(<\/span>n<\/span>):<\/span>\n        # zero down the column<\/span>\n        u<\/span>,<\/span> beta_u<\/span> =<\/span> make_house_vec<\/span>(<\/span>A<\/span>[<\/span>col<\/span>:,<\/span>col<\/span>])<\/span>\n        apply_house_left<\/span>(<\/span>A<\/span>[<\/span>col<\/span>:,<\/span> col<\/span>:],<\/span> u<\/span>,<\/span> beta_u<\/span>,<\/span> U<\/span>,<\/span> m<\/span>)<\/span>\n        if<\/span> col<\/span> <=<\/span> n<\/span>-<\/span>2<\/span>:<\/span>\n            # zero across the row<\/span>\n            v<\/span>,<\/span> beta_v<\/span> =<\/span> make_house_vec<\/span>(<\/span>A<\/span>[<\/span>col<\/span>,<\/span>col<\/span>+<\/span>1<\/span>:]<\/span>.<\/span>T<\/span>)<\/span>\n            apply_house_right<\/span>(<\/span>A<\/span>[<\/span>col<\/span>:,<\/span> col<\/span>+<\/span>1<\/span>:],<\/span> v<\/span>,<\/span> beta_v<\/span>,<\/span> Vt<\/span>,<\/span> n<\/span>)<\/span>\n    return<\/span> U<\/span>,<\/span> A<\/span>,<\/span> Vt<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [5]:<\/div>\n
    \n
    \n
    \n
    A<\/span> =<\/span> np<\/span>.<\/span>arange<\/span>(<\/span>1<\/span>,<\/span>13.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>4<\/span>,<\/span>3<\/span>)<\/span>\nA_test<\/span> =<\/span> np<\/span>.<\/span>copy<\/span>(<\/span>A<\/span>)<\/span>\nprint<\/span> A<\/span>\n\nU<\/span>,<\/span> B<\/span>,<\/span> Vt<\/span> =<\/span> hbd_simple<\/span>(<\/span>A<\/span>)<\/span>\nprint<\/span> "B:<\/span>\\n<\/span>"<\/span>,<\/span> B<\/span>\nprint<\/span> "U:<\/span>\\n<\/span>"<\/span>,<\/span> U<\/span>\nprint<\/span> "Vt:<\/span>\\n<\/span>"<\/span>,<\/span> Vt<\/span>\n\nprint<\/span> "should equal input:"<\/span>,<\/span> np<\/span>.<\/span>allclose<\/span>(<\/span>U<\/span>.<\/span>dot<\/span>(<\/span>B<\/span>)<\/span>.<\/span>dot<\/span>(<\/span>Vt<\/span>),<\/span> A_test<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[  1.   2.   3.]\n [  4.   5.   6.]\n [  7.   8.   9.]\n [ 10.  11.  12.]]\nB:\n[[ 12.88409873  21.87643283   0.        ]\n [  0.           2.24623524  -0.61328133]\n [ -0.           0.           0.        ]\n [ -0.           0.           0.        ]]\nU:\n[[ 0.07761505  0.83305216 -0.54737648  0.01946767]\n [ 0.31046021  0.45123659  0.74434565  0.38203344]\n [ 0.54330537  0.06942101  0.15343813 -0.8224699 ]\n [ 0.77615053 -0.31239456 -0.35040731  0.42096879]]\nVt:\n[[ 1.          0.          0.        ]\n [ 0.          0.66700225  0.7450557 ]\n [ 0.          0.7450557  -0.66700225]]\nshould equal input: True\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Some Optimization<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Now that we have the basic process down, let’s see if we can save some resources. Our major insight is going to be that our input matrix A<\/code>, which ends up with lots of zeros, has lots of room to hold other information. And, to store our rotations, we really only need to store the essential<\/em> part of the Householder vectors. We can do that in the zeroed out areas above and below the diagonal and superdiagonal band. Here’s what it looks like visually:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    \\[
    \n\\begin{bmatrix}
    \nB  &   & \\huge{0} \\\\
    \n   & B &          \\\\
    \n\\huge{0} &  & B   \\\\
    \n & \\huge{0} &
    \n\\end{bmatrix}
    \n\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    And here, we lay out the use of the zero spots for our \\(\\vec{v_r}\\)<\/span> (row vectors) and \\(\\downarrow u_c\\)<\/span>‘s (column vectors). Yes, I think the column vector\/down arrow is very cute. The value after the colon is the length of the vector.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    \\[
    \n\\begin{bmatrix}
    \nd_1 & s_1 & \\vec{v_1}:n-2 & & \\\\
    \n\\downarrow u_1:{m-1} & d_2 & s_2 & \\vec{v_2}:n-3 & \\\\
    \n & \\downarrow u_2:{m-2} & d_2 & s_3 & \\vec{v_3}:n-4 \\\\
    \n & & & \\dots & & & \\\\
    \n & & \\downarrow u_{n-3}:m-(n-3) & d_{n-2} & s_{n-2} & \\vec{v_{n-2}}:1 \\\\
    \n & & & \\downarrow u_{n-2}:m-(n-2)=2  & d_{n-1} & s_{n-1} \\\\
    \n & & & & \\downarrow u_{n-1}:m-(n-1)=1 & d_{n}
    \n\\end{bmatrix}
    \n\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    When we store the information like that, we need a way to extract it. It’s not too bad. If you followed the explict U\/V code above, it’s very similar. We’re simply extracting the \\(U\\)<\/span> and \\(V^T\\)<\/span> back out from the original Householder vectors. There’s one detail with the storage method: we don’t explicitly store the \\(1\\)<\/span> value at the head of the Householder vector.<\/p>\n
    We usually don’t have to explicitly construct \\(U\\)<\/span> or \\(V^T\\)<\/span> to apply them. That is, we can compute the product of \\(U\\)<\/span> or \\(V\\)<\/span> with other matrices efficiently from the packed form. But, we’ll make them explicit here so we can debug our code and check out output.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Extraction<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [6]:<\/div>\n
    \n
    \n
    \n
    # Based on G&VL 3rd. pg. 213 ... 5.1.5<\/span>\ndef<\/span> extractU<\/span>(<\/span>A<\/span>,<\/span> betas<\/span>):<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>\n    Q<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>m<\/span>)<\/span>\n    for<\/span> j<\/span> in<\/span> reversed<\/span>(<\/span>xrange<\/span>(<\/span>n<\/span>)):<\/span> # ignore anything beyond n<\/span>\n        v<\/span> =<\/span> A<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>]<\/span>.<\/span>copy<\/span>()<\/span>\n        v<\/span>[<\/span>0<\/span>]<\/span> =<\/span> 1.0<\/span>\n        Q<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>:]<\/span> =<\/span> (<\/span>np<\/span>.<\/span>eye<\/span>(<\/span>m<\/span>-<\/span>j<\/span>)<\/span> -<\/span> betas<\/span>[<\/span>j<\/span>]<\/span> *<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>v<\/span>,<\/span>v<\/span>))<\/span>.<\/span>dot<\/span>(<\/span>Q<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>:])<\/span>\n    return<\/span> Q<\/span>\n\ndef<\/span> extractVt<\/span>(<\/span>A<\/span>,<\/span> betas<\/span>):<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>\n    Q<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>n<\/span>)<\/span>\n    # stupidly confusing to code b\/c j (our index into the Q<\/span>\n    #                                   we are building)<\/span>\n    # is "off by one"<\/span>\n    # from the values we are extracting in Q<\/span>\n    # this has two practical results:<\/span>\n    #   1.  the last row\/col in Q is ignored (it's an identity row\/col)<\/span>\n    #   2.  to convert from input (A) land to Q land we have to <\/span>\n    #       subtract 1 from the input land stuff (it's all shifted)<\/span>\n    for<\/span> j<\/span> in<\/span> reversed<\/span>(<\/span>xrange<\/span>(<\/span>n<\/span>)):<\/span>\n        v<\/span> =<\/span> A<\/span>[<\/span>j<\/span>-<\/span>1<\/span>,<\/span>j<\/span>+<\/span>1<\/span>-<\/span>1<\/span>:]<\/span>.<\/span>copy<\/span>()<\/span>\n        v<\/span>[<\/span>0<\/span>]<\/span> =<\/span> 1.0<\/span>\n        Q<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>:]<\/span> =<\/span> (<\/span>np<\/span>.<\/span>eye<\/span>(<\/span>n<\/span>-<\/span>j<\/span>)<\/span> -<\/span> betas<\/span>[<\/span>j<\/span>-<\/span>1<\/span>]<\/span> *<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>v<\/span>,<\/span>v<\/span>))<\/span>.<\/span>dot<\/span>(<\/span>Q<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>:])<\/span>\n    return<\/span> Q<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    These two pieces of code are very similar. We can generalize them to a single extract that will pull out either the upper (\\(V^T\\)<\/span>) or lower (\\(U\\)<\/span>) parts.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [7]:<\/div>\n
    \n
    \n
    \n
    def<\/span> extract_house_reflection<\/span>(<\/span>_A<\/span>,<\/span> betas<\/span>,<\/span> side<\/span>=<\/span>"lower"<\/span>):<\/span>\n    A<\/span>,<\/span> shift<\/span> =<\/span> (<\/span>_A<\/span>,<\/span> 0<\/span>)<\/span> if<\/span> side<\/span> ==<\/span> "lower"<\/span> else<\/span> (<\/span>_A<\/span>.<\/span>T<\/span>,<\/span>1<\/span>)<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>\n\n    Q<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>m<\/span>)<\/span>\n    for<\/span> j<\/span> in<\/span> reversed<\/span>(<\/span>xrange<\/span>(<\/span>min<\/span>(<\/span>m<\/span>,<\/span>n<\/span>))):<\/span>\n        v<\/span> =<\/span> A<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>-<\/span>shift<\/span>]<\/span>.<\/span>copy<\/span>()<\/span>\n        v<\/span>[<\/span>0<\/span>]<\/span> =<\/span> 1.0<\/span>\n        miniQ<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>m<\/span>-<\/span>j<\/span>)<\/span> -<\/span> betas<\/span>[<\/span>j<\/span>-<\/span>shift<\/span>]<\/span> *<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>v<\/span>,<\/span>v<\/span>)<\/span>\n        Q<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>:]<\/span> =<\/span> (<\/span>miniQ<\/span>)<\/span>.<\/span>dot<\/span>(<\/span>Q<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>:])<\/span>\n    return<\/span> Q<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    And since we will want to extract all the pieces of the packed result, we can wrap that up in one nice call:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [8]:<\/div>\n
    \n
    \n
    \n
    def<\/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  ... assume rows >= cols<\/span>\n       (works for non-square M)'''<\/span>\n    B<\/span> =<\/span> np<\/span>.<\/span>zeros_like<\/span>(<\/span>M<\/span>)<\/span>\n\n    shape<\/span> =<\/span> B<\/span>.<\/span>shape<\/span>[<\/span>1<\/span>]<\/span>\n    step<\/span> =<\/span> shape<\/span> +<\/span> 1<\/span> # # cols + 1<\/span>\n    end<\/span>  =<\/span> shape<\/span>**<\/span>2<\/span>  # top square (shape,shape)<\/span>\n    \n    B<\/span>.<\/span>flat<\/span>[<\/span> :<\/span>end<\/span>:<\/span>step<\/span>]<\/span> =<\/span> np<\/span>.<\/span>diag<\/span>(<\/span>M<\/span>)<\/span>     # diag<\/span>\n    B<\/span>.<\/span>flat<\/span>[<\/span>1<\/span>:<\/span>end<\/span>:<\/span>step<\/span>]<\/span> =<\/span> np<\/span>.<\/span>diag<\/span>(<\/span>M<\/span>,<\/span>+<\/span>1<\/span>)<\/span>  # super diag<\/span>\n    return<\/span> B<\/span>\n\ndef<\/span> extract_packed_house_bidiag<\/span>(<\/span>H<\/span>):<\/span>\n    U<\/span>  =<\/span> extract_house_reflection<\/span>(<\/span>H<\/span>,<\/span> betas_U<\/span>,<\/span> side<\/span>=<\/span>"lower"<\/span>)<\/span>\n    Vt<\/span> =<\/span> extract_house_reflection<\/span>(<\/span>H<\/span>,<\/span> betas_V<\/span>,<\/span> side<\/span>=<\/span>"upper"<\/span>)<\/span>\n    B<\/span>  =<\/span> extract_upper_bidiag<\/span>(<\/span>H<\/span>)<\/span>\n    return<\/span> U<\/span>,<\/span> B<\/span>,<\/span> Vt<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Implicit Storage<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    And now we turn back to G&VL (3rd) Algorithm 5.4.2 with implicit reflections and reusing A<\/code>‘s storage for u<\/code> and v<\/code>.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [9]:<\/div>\n
    \n
    \n
    \n
    def<\/span> house_bidiag<\/span>(<\/span>A<\/span>):<\/span>\n    m<\/span>,<\/span>n<\/span> =<\/span> A<\/span>.<\/span>shape<\/span>\n    assert<\/span> m<\/span> >=<\/span> n<\/span>\n    U<\/span>,<\/span>V<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>m<\/span>),<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>n<\/span>)<\/span>\n    \n    betas_U<\/span> =<\/span> np<\/span>.<\/span>empty<\/span>(<\/span>n<\/span>)<\/span>\n    betas_V<\/span> =<\/span> np<\/span>.<\/span>empty<\/span>(<\/span>n<\/span>-<\/span>1<\/span>)<\/span>\n    for<\/span> j<\/span> in<\/span> xrange<\/span>(<\/span>n<\/span>):<\/span>\n        u<\/span>,<\/span>betas_U<\/span>[<\/span>j<\/span>]<\/span> =<\/span> make_house_vec<\/span>(<\/span>A<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>])<\/span>\n        miniHouse<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>m<\/span>-<\/span>j<\/span>)<\/span> -<\/span> betas_U<\/span>[<\/span>j<\/span>]<\/span> *<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>u<\/span>,<\/span>u<\/span>)<\/span>\n        A<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>:]<\/span> =<\/span> (<\/span>miniHouse<\/span>)<\/span>.<\/span>dot<\/span>(<\/span>A<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>:])<\/span>\n        # A's rows go 0 to m<\/span>\n        #  A[j,j] and A[j,j+1] are part of the bidiagonal result<\/span>\n        #  A[j+1:,j]  (remainder of this COLUMN) stores u<\/span>\n        #  A[j+1:] means indices j+1, j+2, .... m-1 is m-j positions<\/span>\n        A<\/span>[<\/span>j<\/span>+<\/span>1<\/span>:,<\/span>j<\/span>]<\/span> =<\/span> u<\/span>[<\/span>1<\/span>:]<\/span> # [1:m-j+1]<\/span>\n\n        if<\/span> j<\/span> <<\/span> n<\/span>-<\/span>1<\/span>:<\/span>\n            v<\/span>,<\/span>betas_V<\/span>[<\/span>j<\/span>]<\/span> =<\/span> make_house_vec<\/span>(<\/span>A<\/span>[<\/span>j<\/span>,<\/span>j<\/span>+<\/span>1<\/span>:]<\/span>.<\/span>T<\/span>)<\/span>\n            miniHouse<\/span> =<\/span> np<\/span>.<\/span>eye<\/span>(<\/span>n<\/span>-<\/span>(<\/span>j<\/span>+<\/span>1<\/span>))<\/span> -<\/span> betas_V<\/span>[<\/span>j<\/span>]<\/span> *<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>v<\/span>,<\/span>v<\/span>)<\/span>\n            A<\/span>[<\/span>j<\/span>:,<\/span>j<\/span>+<\/span>1<\/span>:]<\/span> =<\/span> A<\/span>[<\/span>j<\/span>:,<\/span> j<\/span>+<\/span>1<\/span>:]<\/span>.<\/span>dot<\/span>(<\/span>miniHouse<\/span>)<\/span>\n            A<\/span>[<\/span>j<\/span> ,<\/span>j<\/span>+<\/span>2<\/span>:]<\/span> =<\/span> v<\/span>[<\/span>1<\/span>:]<\/span> # [1:n-j]<\/span>\n    return<\/span> betas_U<\/span>,<\/span> betas_V<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    In [10]:<\/div>\n
    \n
    \n
    \n
    A<\/span> =<\/span> np<\/span>.<\/span>arange<\/span>(<\/span>1<\/span>,<\/span>13.0<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>4<\/span>,<\/span>3<\/span>)<\/span>\nA_test<\/span> =<\/span> A<\/span>.<\/span>copy<\/span>()<\/span>\nprint<\/span> A<\/span>\n\nbetas_U<\/span>,<\/span> betas_V<\/span> =<\/span> house_bidiag<\/span>(<\/span>A<\/span>)<\/span>\nU<\/span>,<\/span> B<\/span>,<\/span> Vt<\/span> =<\/span> extract_packed_house_bidiag<\/span>(<\/span>A<\/span>)<\/span>\n\nprint<\/span> " U:<\/span>\\n<\/span>"<\/span>,<\/span> U<\/span>\nprint<\/span> "Vt:<\/span>\\n<\/span>"<\/span>,<\/span> Vt<\/span>\nprint<\/span> " B:<\/span>\\n<\/span>"<\/span>,<\/span> B<\/span>\n\nprint<\/span> "should equal A:"<\/span>,<\/span> np<\/span>.<\/span>allclose<\/span>(<\/span> U<\/span>.<\/span>dot<\/span>(<\/span>B<\/span>)<\/span>.<\/span>dot<\/span>(<\/span>Vt<\/span>),<\/span> A_test<\/span>)<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n
    \n
    <\/div>\n
    \n
    [[  1.   2.   3.]\n [  4.   5.   6.]\n [  7.   8.   9.]\n [ 10.  11.  12.]]\n U:\n[[ 0.07761505  0.83305216 -0.54737648  0.01946767]\n [ 0.31046021  0.45123659  0.74434565  0.38203344]\n [ 0.54330537  0.06942101  0.15343813 -0.8224699 ]\n [ 0.77615053 -0.31239456 -0.35040731  0.42096879]]\nVt:\n[[ 1.          0.          0.        ]\n [ 0.          0.66700225  0.7450557 ]\n [ 0.          0.7450557  -0.66700225]]\n B:\n[[ 12.88409873  21.87643283   0.        ]\n [  0.           2.24623524  -0.61328133]\n [  0.           0.           0.        ]\n [  0.           0.           0.        ]]\nshould equal A: True\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":"
    We’re going to make use of Householder reflections to turn a generic matrix into a bidiagonal matrix. I’ve laid a lot of foundation for these topics in the posts I just linked. Check them out! Onward.<\/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-445","post","type-post","status-publish","format-standard","hentry","category-mrdr","category-sci-math-stat-python"],"_links":{"self":[{"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/445","targetHints":{"allow":["GET"]}}],"collection":[{"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/comments?post=445"}],"version-history":[{"count":1,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/445\/revisions"}],"predecessor-version":[{"id":446,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/445\/revisions\/446"}],"wp:attachment":[{"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/media?parent=445"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/categories?post=445"},{"taxonomy":"post_tag","embeddable":true,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/tags?post=445"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}