{"id":456,"date":"2017-02-02T19:35:23","date_gmt":"2017-02-02T19:35:23","guid":{"rendered":"http:\/\/drsfenner.org\/blog\/?p=456"},"modified":"2017-02-02T19:35:40","modified_gmt":"2017-02-02T19:35:40","slug":"dc-svd-i-just-cant-let-it-go","status":"publish","type":"post","link":"http:\/\/drsfenner.org\/blog\/2017\/02\/dc-svd-i-just-cant-let-it-go\/","title":{"rendered":"DC SVD I: Just Can’t Let It Go"},"content":{"rendered":"
\n
\n<\/div>\n
\n
\n

It’s been way, way, way too long since I’ve posted. I haven’t been slacking though, I’ve merely been busy. Really.<\/p>\n

I decided to dive back into the SVD problem and look at an alternative to the QR based SVD computations. Namely, I’m going to give a breakdown of the divide-and-conquer<\/em> approach to computing the SVD. Similar to the relationship between bidiagonal SVD and tridiagonal QR decompositions, there is a close relationship between dividing-and-conquering QR and SVD. I’m going to start from "the inside out" with the innermost task of DC (divide-and-conquer) process: solving a particular equation known as the secular equation (secular meaning "not heavenly" – i.e., earthly or planetly – check it out on wikipedia).<\/p>\n

One note: I feel very guilty about still using Python 2. I had intented to transition this project to Python 3 over the course of the last year. Alas, my major training clients over the last year were using Python 2 and I didn’t really want to have a mixed development environment on my laptop. Well, at least there’s something to do if I make a book out of these posts!<\/p>\n

Enough preamble, let’s get down to business.<\/p>\n<\/div>\n<\/div>\n<\/div>\n

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

The Secular Equation<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n

Divide-and-conquer SVD is built on computing the roots of the secular equation<\/em>. The roots, or zeros, of an equation are the \\(\\lambda\\)<\/span>s such that \\(f(\\lambda) = 0\\)<\/span>. The secular equation is (where we take \\(\\rho=1\\)<\/span>): \\[f(\\lambda) = 1 + \\rho \\sum_{i=1}^n \\frac{z_{i}^2}{d_i – \\lambda} = 1 + \\sum_{i=1}^n \\frac{z_{i}^2}{d_i – \\lambda}\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n

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

Here is a graph of a secular equation and its roots:<\/p>\n<\/div>\n<\/div>\n<\/div>\n

\n
\n
In [1]:<\/div>\n
\n
\n
\n
<\/span>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 \nxs<\/span> =<\/span> np<\/span>.<\/span>linspace<\/span>(<\/span>-<\/span>1<\/span>,<\/span>10<\/span>,<\/span>10000<\/span>)<\/span>\n\n# graph secular equation (blue curve)<\/span>\n# plot z,d --> y<\/span>\nzs<\/span> =<\/span> np<\/span>.<\/span>array<\/span>([<\/span>.<\/span>6<\/span>,<\/span>1.2<\/span>,<\/span>1.8<\/span>])<\/span>\nds<\/span> =<\/span> np<\/span>.<\/span>array<\/span>([<\/span>0.0<\/span>,<\/span>3<\/span>,<\/span>5<\/span>])<\/span>\nwith<\/span> np<\/span>.<\/span>errstate<\/span>(<\/span>divide<\/span>=<\/span>'ignore'<\/span>):<\/span>\n    # plain vanilla to "fancy-schmancy"<\/span>\n    #ys1 = 1.0 + (1.0 \/ (0.0-xs_sq)) + (4.0 \/ (9-xs_sq)) + (16.0 \/ (25 - xs_sq))<\/span>\n    #ys2 = 1.0 + (zs[0]**2 \/ (ds[0]**2-xs_sq)) + (zs[1]**2 \/ (ds[1]**2-xs_sq)) + (zs[2]**2 \/ (ds[2]**2 - xs_sq))<\/span>\n    ys3<\/span> =<\/span> 1.0<\/span> +<\/span> ((<\/span>zs<\/span>**<\/span>2<\/span>)<\/span>.<\/span>reshape<\/span>(<\/span>3<\/span>,<\/span>1<\/span>)<\/span> \/<\/span> np<\/span>.<\/span>subtract<\/span>.<\/span>outer<\/span>(<\/span>ds<\/span>,<\/span> xs<\/span>))<\/span>.<\/span>sum<\/span>(<\/span>0<\/span>)<\/span>\n    # assert np.allclose(ys1, ys2) and np.allclose(ys2,ys3)<\/span>\nplt<\/span>.<\/span>plot<\/span>(<\/span>xs<\/span>,<\/span> ys3<\/span>)<\/span>\n\n# add an x-axis (yellow horizontal line)<\/span>\nplt<\/span>.<\/span>plot<\/span>(<\/span>xs<\/span>,<\/span> np<\/span>.<\/span>zeros_like<\/span>(<\/span>xs<\/span>),<\/span> 'y-'<\/span>)<\/span>\n\n# add poles\/asymptotes (grey vertical lines)<\/span>\nplt<\/span>.<\/span>vlines<\/span>(<\/span>ds<\/span>,<\/span> -<\/span>10<\/span>,<\/span> 10<\/span>,<\/span> '.75'<\/span>)<\/span>\n\n# add roots (red dots)<\/span>\n# use equivalent matrix (for these z,d) and eigenvalue computation to find roots<\/span>\n# ds are diagonal entries, zs are first columns (note, ds[0] is 0.0 by definition)<\/span>\nROU<\/span> =<\/span> np<\/span>.<\/span>diag<\/span>(<\/span>ds<\/span>)<\/span> +<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>zs<\/span>,<\/span> zs<\/span>)<\/span>\n#tm[:,0] = zs<\/span>\nzeros_act<\/span> =<\/span> nla<\/span>.<\/span>eig<\/span>(<\/span>ROU<\/span>)[<\/span>0<\/span>]<\/span>\nplt<\/span>.<\/span>plot<\/span>(<\/span>zeros_act<\/span>,<\/span> np<\/span>.<\/span>zeros_like<\/span>(<\/span>zeros_act<\/span>),<\/span> 'r.'<\/span>)<\/span>\n\n#scaling<\/span>\nplt<\/span>.<\/span>ylim<\/span>(<\/span>-<\/span>10<\/span>,<\/span>10<\/span>),<\/span> plt<\/span>.<\/span>xlim<\/span>(<\/span>-<\/span>1<\/span>,<\/span> 10<\/span>);<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n
\n
<\/div>\n
\n\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
A quick note, the secular equation as written above is most directly used to compute eigenvalues<\/em> not singular values<\/em>. There is a strong relationship between the two sets of values and they lead to only slightly different forms in the secular equation. I mention this because you might see slightly different forms of the secular equation depending on whether you are reading about eigenvalues or singular values. We will use the form above for solving both problems by slightly modifying the resulting \\(\\lambda\\)<\/span>s.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Back to our regularly scheduled program. We will find our desired roots by isolating out \\(f(\\lambda)\\)<\/span> between each set of poles. The poles occur at the values of \\(d_i\\)<\/span>. So, on the interval \\((d_i, d_{i+1})\\)<\/span>, the problem simplifies to finding a single root of the secular equation. Throughout this post, we’ll assume that \\(d_i < d_j\\)<\/span> for \\(i<j\\)<\/span> (in English, the \\(d_i\\)<\/span>s are distinct and sorted). We’ll find a single root \\(n\\)<\/span> times and find all \\(n\\)<\/span> roots. The cleverest will now point out that \\(n\\)<\/span> poles only hold \\(n-1\\)<\/span> zeros between them. The last zero is to the right of \\(d_n\\)<\/span>.<\/p>\n
How do we find the single roots?<\/p>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Newton’s Method<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
\n
\n<\/div>\n
\n
\n
Let’s take a second and review Newton’s method<\/em> for finding a root of an equation \\(f(x)\\)<\/span> near \\(x_0\\)<\/span>.<\/p>\n
    \n
  1. Approximate \\(f(x)\\)<\/span> by a linear function \\(l(x) = ax+b\\)<\/span><\/li>\n
  2. Apply constraints such that \\(f(x_0) = l(x_0)\\)<\/span> and \\(f'(x_0) = l'(x_0)\\)<\/span>. These determine the coefficients of \\(l(x)\\)<\/span>.<\/li>\n
  3. Find the root of \\(l(x)\\)<\/span> and call it \\(x_1\\)<\/span>. In this case, the root is the \\(x\\)<\/span>-intercept of \\(l(x)\\)<\/span>.<\/li>\n
  4. That root of \\(l(x)\\)<\/span> is a better guess as to the root of \\(f(x)\\)<\/span>.<\/li>\n<\/ol>\n
    Newton’s method is a great technique and it is used broadly because it is conceptually and formally simple, easy to implement, and the iteration steps are reasonably fast. However, in the case of the secular equation, as the numerators get very small, the corner in the graph will go from a gentle bend to a sharp corner. In other words, the graph will go from vertical to horizontal very quickly. When this happens, a linear approximation on the nearly horizontal part (which is the majority of the pole-to-pole interval), will also be approximately horizontal and be aimed far away from the true root in this interval.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    A Modified Newton’s Method<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    So, we need to try something else. Fortunately, we can maintain the outline of Newton’s method while using a different approximating function. Instead of using a linear form, we will use the following rational function of \\(\\lambda\\)<\/span> which has poles at \\(d_i\\)<\/span> and \\(d_{i+1}\\)<\/span>:<\/p>\n
    \\[h(\\lambda) = \\frac{C_1}{d_i-\\lambda} + \\frac{C_2}{d_{i+1}-\\lambda} + C_3\\]<\/span><\/p>\n
    So, we are approximating \\(f(\\lambda)\\)<\/span> with \\(h(\\lambda)\\)<\/span>:<\/p>\n
    \\[f(\\lambda) = 1 + \\sum_{i=1}^n \\frac{z_{i}^2}{d_i – \\lambda} \\approx
    \n\\frac{C_1}{d_i-\\lambda} + \\frac{C_2}{d_{i+1}-\\lambda} + C_3 = h(x)\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Also, since we want to avoid numerical problems with term cancellation, we will break the sum in \\(f(\\lambda)\\)<\/span> into two parts: (1) the sum up to term \\(k\\)<\/span> and (2) the sum from term \\(k+1\\)<\/span> on. Along with the an ordering assumption on the \\(d_i\\)<\/span>, this means that the lower terms are all negative and the upper terms are all positive for \\(\\lambda \\in (d_i, d_{i+1})\\)<\/span>. We will also give names to those partial sums.<\/p>\n
    \\[
    \n\\begin{eqnarray}
    \nf(\\lambda) &=& 1 + & \\sum_{i=1}^k \\frac{z_{i}^2}{d_i – \\lambda} + \\sum_{i=k+1}^n \\frac{z_{i}^2}{d_i – \\lambda} &=& 1 + \\Psi_1(\\lambda) + \\Psi_2(\\lambda)\\\\
    \nf'(\\lambda) &=& & \\sum_{i=1}^k \\frac{z_i^2}{(d_i – \\lambda)^2} + \\sum_{i=k+1}^n \\frac{z_i^2}{(d_i – \\lambda)^2} &=& \\Psi'_1(\\lambda) + \\Psi'_2(\\lambda)
    \n\\end{eqnarray}
    \n\\]<\/span><\/p>\n
    Note that the derivative (and the derivatives of the partial sums) is strictly positive except at \\(\\lambda = d_i\\)<\/span> (the poles of \\(f\\)<\/span>).<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Since we have broken \\(f\\)<\/span> into pieces, we will also break \\(h\\)<\/span> into pieces:<\/p>\n
    \\[h(\\lambda)=1 + h_1(\\lambda) + h_2(\\lambda)\\]<\/span><\/p>\n
    and we will give them each the same form, which is similar to that of \\(h\\)<\/span>:<\/p>\n
    \\[h_1(\\lambda) = \\hat{c}_1 + \\frac{C_1}{d_k – \\lambda} \\quad
    \nh_2(\\lambda) = \\hat{c}_2 + \\frac{C_2}{d_{k+1} – \\lambda}\\]<\/span><\/p>\n
    We will also enforce the "Newton conditions"<\/em> that both the value of the approximations \\(h_i\\)<\/span> and the value of the derivative of the approximations \\(h'_i\\)<\/span> agree with the target function at \\(\\lambda_{\\alpha}\\)<\/span>:<\/p>\n
    \\[
    \n\\begin{align}
    \nh_1(\\lambda_{\\alpha}) = \\Psi_1(\\lambda_{\\alpha}) &\\quad& h'_1(\\lambda_{\\alpha}) = \\Psi'_1(\\lambda_{\\alpha})\\\\
    \nh_2(\\lambda_{\\alpha}) = \\Psi_2(\\lambda_{\\alpha}) &\\quad& h'_2(\\lambda_{\\alpha}) = \\Psi'_2(\\lambda_{\\alpha})
    \n\\end{align}
    \n\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    The result of these constraints is that: \\[ C_1 = \\Psi'_1(\\lambda_{\\alpha})(d_k – \\lambda_{\\alpha})^2 \\\\
    \n   \\hat{c}_1 = \\Psi_1(\\lambda_{\\alpha}) – \\Psi'_1(\\lambda_{\\alpha})(d_k – \\lambda_{\\alpha})\\]<\/span><\/p>\n
    \\[ C_2 = \\Psi'_2(\\lambda_{\\alpha})(d_k – \\lambda_{\\alpha})^2 \\\\
    \n   \\hat{c}_2 = \\Psi_2(\\lambda_{\\alpha}) – \\Psi'_2(\\lambda_{\\alpha})(d_{k+1} – \\lambda_{\\alpha})\\]<\/span><\/p>\n
    \\[C_3 = 1 + \\hat{c}_1 + \\hat{c}_2\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Since we are seeking a root of \\(h(\\lambda)\\)<\/span>, we can set<\/p>\n
    \\[h(\\lambda) = \\frac{C_1}{d_i-\\lambda} + \\frac{C_2}{d_{i+1}-\\lambda} + C_3 = 0\\]<\/span><\/p>\n
    and multiply out the denominators. This gives us:<\/p>\n
    \\[0 = C_1 (d_{k+1} – \\lambda) + C_2 (d_k – \\lambda) + C_3(d_k – \\lambda)(d_{k+1} – \\lambda) \\qquad \\star\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    In our original Newton formulation of root finding, we went from a guess of the root \\(x_0\\)<\/span> to \\(x_1\\)<\/span>. We could phrase that as going \\(x_1=x_0+\\kappa\\)<\/span> where \\(\\kappa\\)<\/span> is a correction to our our \\(x_1\\)<\/span> guess. We will use that "\\(\\kappa\\)<\/span>-technique" here. Call our current approximation or guess \\(\\lambda_\\alpha\\)<\/span>. Let’s call our correction \\(\\eta\\)<\/span>. Then we can write a \\(\\lambda_{\\alpha} + \\eta = \\lambda\\)<\/span>. Another way to think about this is as a variable substitution (or function shift) for \\(\\lambda \\rightarrow \\lambda – \\lambda_{\\alpha} = \\eta\\)<\/span> which centers the function horizontally at \\(\\lambda_{\\alpha}\\)<\/span>. If we solve for a \\(\\lambda\\)<\/span> that gives a zero, we would also have a corresponding \\(\\eta\\)<\/span> value where \\(\\lambda_{\\alpha} + \\eta\\)<\/span> gives a zero. Also, the equality means we can write \\(\\lambda_\\alpha = \\lambda – \\eta\\)<\/span>.<\/p>\n
    Now, to help simplify the terms, let \\(\\Delta_k = d_k – \\lambda_\\alpha\\)<\/span> and \\(\\Delta_{k+1} = d_{k+1} – \\lambda_\\alpha\\)<\/span>. Together, these give us (and likewise for \\(\\Delta_{k+1}\\)<\/span>):<\/p>\n
    \\[\\Delta_k = d_k – \\lambda_\\alpha = d_k – (\\lambda – \\eta) = d_k – \\lambda + \\eta\\]<\/span><\/p>\n
    Rearranging:<\/p>\n
    \\[d_k – \\lambda = \\Delta_k – \\eta \\quad \\text{and} \\quad d_{k+1} – \\lambda = \\Delta_{k+1} – \\eta\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Since, \\(\\lambda_\\alpha\\)<\/span> is a constant during an interation of Newton’s method, we will drop the functional depence on it. Our \\(C\\)<\/span> constants become: \\[\\begin{align}
    \nC_1 &= \\Psi'_1 \\Delta_{k}^2 \\\\
    \nC_2 &= \\Psi'_2 \\Delta_{k+1}^2 \\\\
    \nC_3 &= 1 + \\Psi_1 + \\Psi_2 – \\Psi'_1 \\Delta_{k} – \\Psi'_2 \\Delta_{k+1}
    \n\\end{align}
    \n\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    The punchline is an expression we can solve for \\(\\eta\\)<\/span> (the correction to our old guess \\(\\lambda_\\alpha\\)<\/span>):<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    \\[\\begin{align}
    \nC_1 (d_{k+1} – \\lambda) + C_2 (d_k – \\lambda) + C_3(d_k – \\lambda)(d_{k+1} – \\lambda)
    \n    &= C_1 (\\Delta_{k+1} – \\eta) + C_2 (\\Delta_{k} – \\eta) + C_3 (\\Delta_{k} – \\eta) (\\Delta_{k+1} – \\eta) \\\\
    \n    &= {\\eta}^2 C_3\\ + \\\\
    \n    &\\quad {-\\eta} (C_1 + C_2 + C_3(\\Delta_k + \\Delta_{k+1}))\\ +\\\\
    \n    &{\\quad}C_1\\Delta_{k+1} + C_2\\Delta_{k} + C_3\\Delta_k\\Delta_{k+1}
    \n\\end{align}
    \n\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    So, rewriting these constants into the form of the classic quadratic equation \\(a\\eta^2+b\\eta+c\\)<\/span>, we have (we introduce \\(\\Psi_{\\text{sum}} = 1 + \\Psi_1 + \\Psi_2\\)<\/span>):<\/p>\n
    \\[
    \n\\begin{align}
    \na &= C_3 \\\\
    \n  &= (1 + \\Psi_1 + \\Psi_2) – \\Psi'_1 \\Delta_{k} – \\Psi'_2 \\Delta_{k+1} \\\\
    \n  &= \\Psi_{\\text{sum}} – \\Psi'_1 \\Delta_{k} – \\Psi'_2 \\Delta_{k+1} \\\\
    \n\\\\
    \n-b &= C_1 + C_2 + C_3(\\Delta_k + \\Delta_{k+1}) \\\\
    \n&= \\Psi'_1 \\Delta_{k}^2 + \\Psi'_2 \\Delta_{k+1}^2 + C_3(\\Delta_k + \\Delta_{k+1}) \\\\
    \n&= \\Psi'_1 \\Delta_{k}^2 + \\Psi'_2 \\Delta_{k+1}^2 + (\\Psi_{\\text{sum}} – \\Psi'_1 \\Delta_{k} – \\Psi'_2 \\Delta_{k+1})(\\Delta_k + \\Delta_{k+1})\\\\
    \n&= \\Psi_{\\text{sum}} (\\Delta_k + \\Delta_{k+1}) – \\Psi'_1 \\Delta_{k} \\Delta_{k+1} – \\Psi'_1 \\Delta_{k} \\Delta_{k+1} \\\\
    \n&= \\Psi_{\\text{sum}} (\\Delta_k + \\Delta_{k+1}) – (\\Psi'_1 + \\Psi'_2)\\Delta_{k} \\Delta_{k+1} \\\\
    \nb &= (\\Psi'_1 + \\Psi'_2)\\Delta_{k} \\Delta_{k+1} – \\Psi_{\\text{sum}} (\\Delta_k + \\Delta_{k+1}) \\\\
    \n\\\\
    \nc &= C_3\\Delta_k\\Delta_{k+1} + C_1\\Delta_{k+1} + C_2\\Delta_{k} \\\\
    \n&= (\\Psi_{\\text{sum}} – \\Psi'_1 \\Delta_{k} – \\Psi'_2 \\Delta_{k+1})\\Delta_k\\Delta_{k+1} + \\Psi'_1 \\Delta_{k}^2\\Delta_{k+1} + \\Psi'_2 \\Delta_{k+1}^2 \\Delta_{k} \\\\
    \n&= \\Psi_{\\text{sum}} \\Delta_{k} \\Delta_{k+1}
    \n\\end{align}
    \n\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    The solutions \\(\\eta_1, \\eta_2\\)<\/span> to the quadratic above are well known as \\(\\eta_1, \\eta_2 = \\frac{-b \\pm \\sqrt{b^2-4ac}}{2a} = \\frac{-b \\pm \\text{root}}{2a}\\)<\/span>. However, there is an interesting (and useful!) alternative formula for the solutions: \\(\\eta_2,\\eta_1 = \\frac{2c}{-b \\mp \\sqrt{b^2-4ac}} = \\frac{2c}{-b \\mp \\text{root}}\\)<\/span>. These come from either (1) using different steps to solve the generic quadratic form or (2) applying Vieta’s formulas to the standard solution. Due to some algebraic magic and inequalities I haven’t quite figured out yet, the solution that we need can be picked out based on the value of \\(b\\)<\/span>. (Li claims – and tests on running code show – that based on the quadratic equation and the constaints above – including \\(C_1,C_2 > 0\\)<\/span> – the correct \\(\\eta\\)<\/span> is based on \\(b\\)<\/span>. I can’t demonstrate that, despite trying many different algebraic contortions. I did get to the point of showing \\(C_3>0\\)<\/span>, but alas. I’d love help and I’ll buy you a cookie or a beer for it.) The punchline is:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    \\[
    \n\\eta=\\begin{cases}
    \n\\frac{2c}{-b + \\text{root}} & b<0 \\\\
    \n\\frac{-b – \\text{root}}{2a} & b \\geq 0
    \n\\end{cases}
    \n\\]<\/span><\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Two last pieces lead us to our improvised (MacGyvered?) Newton’s method.<\/p>\n
    We need to know where to start and when to stop. These topics make varied appearances in the literature. Since all but one root falls between the poles, we can start most root finding iterations at the midpoint of the poles \\(\\frac{d_i + d_{i+1}}{2}\\)<\/span>. For the straggler root (which lies to the right of \\(d_n\\)<\/span>), we use a starting point at \\(d_n + \\|z\\|^2\\)<\/span> (following Arbenz’s lead). I don’t use a midpoint to the right (i.e., \\(\\frac{d_n + (d_n + \\|z\\|^2)}{2} = d_n + \\frac{\\|z\\|^2}{2}\\)<\/span>) because I’ve had small numerical differences with LAPACK’s results. Another note, Gu & Eisenstat say the root should be to the left of \\(d_n + \\|z\\|\\)<\/span>. So, there’s some discreprancy here. Oh well. It happens.<\/p>\n
    Gu & Eisenstat<\/a> propose stopping when \\(|\\Psi_{\\text{sum}}| < 10.0\\epsilon(1 + |\\Psi_1| + |\\Psi_2|)\\)<\/span>.<\/p>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    Code Outline<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
    \n
    \n<\/div>\n
    \n
    \n
    For all that ridiculously long development, the code becomes:<\/p>\n
      \n
    1. Setup: set some variables (including inital \\(\\lambda\\)<\/span>), deal with \\(n\\)<\/span>-th root, determine closest \\(d_i\\)<\/span>, and shift (center) the problem there<\/li>\n
    2. Loop until stopping condition is met\n
        \n
      1. Calculate \\(\\Psi\\)<\/span>s, \\(\\Delta\\)<\/span>s, \\(a, b, c\\)<\/span>, and solve the quadratic for \\(\\eta\\)<\/span>.<\/li>\n
      2. Update \\(\\lambda\\)<\/span>.<\/li>\n<\/ol>\n<\/li>\n
      3. Undo the shift.<\/li>\n
      4. Et voila.<\/li>\n<\/ol>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        In [2]:<\/div>\n
        \n
        \n
        \n
        <\/span>import<\/span> itertools<\/span> as<\/span> it<\/span>\n\nEPSILON<\/span> =<\/span> np<\/span>.<\/span>finfo<\/span>(<\/span>float<\/span>)<\/span>.<\/span>eps<\/span>\nNEPS<\/span> =<\/span> 10.0<\/span> *<\/span> EPSILON<\/span> # via Gu & Eisenstat<\/span>\n\ndef<\/span> solve_secular_oneroot<\/span>(<\/span>d<\/span>,<\/span> z<\/span>,<\/span> i<\/span>):<\/span>\n    # translated and modified from:  Arbenz, Peter. see citation below.<\/span>\n    # there are no "l" (ells) in the following code, they are all "1" (ones)<\/span>\n    n<\/span> =<\/span> d<\/span>.<\/span>size<\/span>;<\/span> d_i<\/span> =<\/span> d<\/span>[<\/span>i<\/span>];<\/span> z<\/span> =<\/span> z<\/span>*<\/span>z<\/span>\n    \n    if<\/span> i<\/span> <<\/span> n<\/span> -<\/span> 1<\/span>:<\/span>\n        d_ip1<\/span>   =<\/span> d<\/span>[<\/span>i<\/span>+<\/span>1<\/span>]<\/span>\n        lambda_<\/span> =<\/span> (<\/span>d_i<\/span> +<\/span> d_ip1<\/span>)<\/span> \/<\/span> 2.0<\/span>\n    else<\/span>:<\/span>\n        adder<\/span> =<\/span> max<\/span>(<\/span>nla<\/span>.<\/span>norm<\/span>(<\/span>z<\/span>)<\/span>**<\/span>2<\/span>,<\/span> 1.0<\/span>)<\/span> # protect against small norm<\/span>\n        d_ip1<\/span>   =<\/span> d<\/span>[<\/span>-<\/span>1<\/span>]<\/span> +<\/span> adder<\/span>\n        lambda_<\/span> =<\/span> d_ip1<\/span>\n\n    left_idx<\/span>,<\/span> right_idx<\/span> =<\/span> slice<\/span>(<\/span>0<\/span>,<\/span>i<\/span>+<\/span>1<\/span>),<\/span> slice<\/span>(<\/span>i<\/span>+<\/span>1<\/span>,<\/span> None<\/span>)<\/span>\n    left_z<\/span>,<\/span> right_z<\/span> =<\/span> z<\/span>[<\/span>left_idx<\/span>],<\/span> z<\/span>[<\/span>right_idx<\/span>]<\/span>\n\n    denoms<\/span> =<\/span> d<\/span> -<\/span> lambda_<\/span>    # PN: this makes a new vector<\/span>\n    psi1<\/span> =<\/span> np<\/span>.<\/span>sum<\/span>(<\/span> left_z<\/span>**<\/span>2<\/span> \/<\/span> denoms<\/span>[<\/span> left_idx<\/span>])<\/span>\n    psi2<\/span> =<\/span> np<\/span>.<\/span>sum<\/span>(<\/span>right_z<\/span>**<\/span>2<\/span> \/<\/span> denoms<\/span>[<\/span>right_idx<\/span>])<\/span>\n\n    # determine whether soln is closer to right or left pole<\/span>\n    # (remember, we are currently at lambda == midpt)<\/span>\n    psi_sum<\/span> =<\/span> 1.0<\/span> +<\/span> psi1<\/span> +<\/span> psi2<\/span>\n    left_half<\/span> =<\/span> psi_sum<\/span> ><\/span> 0.0<\/span>\n    active_d<\/span> =<\/span> d_i<\/span> if<\/span> left_half<\/span> else<\/span> d_ip1<\/span>\n\n    # shift poles (so one of them is zero), initial guess, and d_i,d_ip1<\/span>\n    d<\/span> =<\/span> d<\/span> -<\/span> active_d<\/span>\n    lambda_<\/span> =<\/span> lambda_<\/span> -<\/span> active_d<\/span>\n\n    if<\/span> left_half<\/span>:<\/span>\n        d_i<\/span>,<\/span> d_ip1<\/span> =<\/span> 0.0<\/span>,<\/span> d_ip1<\/span>-<\/span>d_i<\/span>\n    else<\/span>:<\/span>\n        d_i<\/span>,<\/span> d_ip1<\/span> =<\/span> d_i<\/span>-<\/span>d_ip1<\/span>,<\/span> 0.0<\/span>\n\n    stop_cond<\/span> =<\/span> NEPS<\/span> *<\/span> n<\/span> *<\/span> (<\/span>1.0<\/span> +<\/span> np<\/span>.<\/span>abs<\/span>(<\/span>psi1<\/span>)<\/span> +<\/span> np<\/span>.<\/span>abs<\/span>(<\/span>psi2<\/span>))<\/span>\n    itct<\/span> =<\/span> it<\/span>.<\/span>count<\/span>()<\/span>\n    while<\/span> (<\/span>np<\/span>.<\/span>abs<\/span>(<\/span>psi_sum<\/span>)<\/span> ><\/span> stop_cond<\/span> and<\/span> itct<\/span>.<\/span>next<\/span>()<\/span> <<\/span> 10<\/span>):<\/span>\n        denoms<\/span> =<\/span> d<\/span> -<\/span> lambda_<\/span>\n        \n        psi1<\/span>,<\/span> psi1_p<\/span> =<\/span> (<\/span>np<\/span>.<\/span>sum<\/span>(<\/span>left_z<\/span> \/<\/span> denoms<\/span>[<\/span> left_idx<\/span>]),<\/span>\n                        np<\/span>.<\/span>sum<\/span>(<\/span>left_z<\/span> \/<\/span> denoms<\/span>[<\/span> left_idx<\/span>]<\/span>**<\/span>2<\/span>))<\/span>\n        psi2<\/span>,<\/span> psi2_p<\/span> =<\/span> (<\/span>np<\/span>.<\/span>sum<\/span>(<\/span>right_z<\/span> \/<\/span> denoms<\/span>[<\/span>right_idx<\/span>]),<\/span>\n                        np<\/span>.<\/span>sum<\/span>(<\/span>right_z<\/span> \/<\/span> denoms<\/span>[<\/span>right_idx<\/span>]<\/span>**<\/span>2<\/span>))<\/span>\n        psi_sum<\/span> =<\/span> 1.0<\/span> +<\/span> psi1<\/span> +<\/span> psi2<\/span>\n\n\n        # D is \\Delta : recall, we've shifted d_i, d_ip1 (line 36 above)<\/span>\n        if<\/span> left_half<\/span>:<\/span>\n            D_i<\/span>,<\/span> D_ip1<\/span> =<\/span> -<\/span>lambda_<\/span>,<\/span> d_ip1<\/span> -<\/span> lambda_<\/span>\n        else<\/span>:<\/span>\n            D_i<\/span>,<\/span> D_ip1<\/span> =<\/span> d_i<\/span> -<\/span> lambda_<\/span>,<\/span> -<\/span>lambda_<\/span>\n\n        #a \\eta^2 + b \\eta + c = 0<\/span>\n        a<\/span> =<\/span> psi_sum<\/span> -<\/span> (<\/span>D_i<\/span> *<\/span> psi1_p<\/span>)<\/span> -<\/span> (<\/span>D_ip1<\/span> *<\/span> psi2_p<\/span>)<\/span>\n        b<\/span> =<\/span> -<\/span> (<\/span>D_i<\/span> +<\/span> D_ip1<\/span>)<\/span> *<\/span> psi_sum<\/span> +<\/span> (<\/span>D_i<\/span> *<\/span> D_ip1<\/span> *<\/span> (<\/span>psi1_p<\/span> +<\/span> psi2_p<\/span>))<\/span>\n        c<\/span> =<\/span> D_i<\/span> *<\/span> D_ip1<\/span> *<\/span> psi_sum<\/span>\n\n        root<\/span> =<\/span> np<\/span>.<\/span>sqrt<\/span>(<\/span>b<\/span>**<\/span>2<\/span> -<\/span> 4<\/span>*<\/span>a<\/span>*<\/span>c<\/span>)<\/span>\n        eta<\/span> =<\/span> (<\/span>2.0<\/span> *<\/span> c<\/span>)<\/span> \/<\/span> (<\/span>-<\/span>b<\/span> +<\/span> root<\/span>)<\/span> if<\/span> b<\/span> <<\/span> 0.0<\/span> else<\/span> (<\/span>-<\/span>b<\/span> -<\/span> root<\/span>)<\/span>\/<\/span>(<\/span>2.0<\/span>*<\/span>a<\/span>)<\/span>\n        lambda_<\/span> =<\/span> lambda_<\/span> +<\/span> eta<\/span>\n\n        stop_cond<\/span> =<\/span> NEPS<\/span> *<\/span> n<\/span> *<\/span> (<\/span>1.0<\/span> +<\/span> np<\/span>.<\/span>abs<\/span>(<\/span>psi1<\/span>)<\/span> +<\/span> np<\/span>.<\/span>abs<\/span>(<\/span>psi2<\/span>))<\/span>\n\n    d_m_lambda<\/span> =<\/span> d<\/span> -<\/span> lambda_<\/span>      # (both shifted) - this will be discussed next post<\/span>\n    lambda_<\/span> =<\/span> lambda_<\/span> +<\/span> active_d<\/span>\n    return<\/span> lambda_<\/span>,<\/span> d_m_lambda<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        Testing<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        I’d be remiss, if I didn’t offer some code to test the root-finding code. The basic process is to use NumPy’s (and hence LAPACK’s) eigenvalue computation system to calculate the roots in a "known good" fashion. NumPy wants a matrix as input, so we take our \\(d\\)<\/span> and \\(z\\)<\/span> values and put them in an appropriate matrix. The actual test cases here are a bit ad hoc<\/em> – they came up during development of the larger D&C SVD code. I’ll have more thorough testing either in the next post or as a standalone follow-up to the next post.<\/p>\n
        Here goes:<\/p>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        In [4]:<\/div>\n
        \n
        \n
        \n
        <\/span>np<\/span>.<\/span>set_printoptions<\/span>(<\/span>precision<\/span>=<\/span>10<\/span>,<\/span>suppress<\/span>=<\/span>False<\/span>)<\/span>\n\ndef<\/span> run_for<\/span>(<\/span>d<\/span>,<\/span>z<\/span>,<\/span> verbose<\/span>=<\/span>False<\/span>):<\/span>\n    assert<\/span> d<\/span>[<\/span>0<\/span>]<\/span> ==<\/span> 0.0<\/span> and<\/span> d<\/span>.<\/span>shape<\/span> ==<\/span> z<\/span>.<\/span>shape<\/span>\n    # ROU is a "rank-one updated" matrix<\/span>\n    ROU<\/span> =<\/span> np<\/span>.<\/span>diag<\/span>(<\/span>d<\/span>)<\/span> +<\/span> np<\/span>.<\/span>outer<\/span>(<\/span>z<\/span>,<\/span>z<\/span>)<\/span>\n    if<\/span> verbose<\/span>:<\/span>\n        print<\/span> "cond(ROU): <\/span>%5.2g<\/span>"<\/span> %<\/span> nla<\/span>.<\/span>cond<\/span>(<\/span>ROU<\/span>)<\/span>\n        print<\/span> "Effective ROU:<\/span>\\n<\/span>"<\/span>,<\/span> ROU<\/span>\n    expected<\/span> =<\/span> np<\/span>.<\/span>sort<\/span>(<\/span>nla<\/span>.<\/span>eigvals<\/span>(<\/span>ROU<\/span>))<\/span>\n    sols<\/span> =<\/span> [<\/span>solve_secular_oneroot<\/span>(<\/span>d<\/span>,<\/span>z<\/span>,<\/span>i<\/span>)[<\/span>0<\/span>]<\/span> for<\/span> i<\/span> in<\/span> xrange<\/span>(<\/span>d<\/span>.<\/span>size<\/span>)]<\/span>\n    actual<\/span> =<\/span> np<\/span>.<\/span>sort<\/span>(<\/span>np<\/span>.<\/span>array<\/span>(<\/span>sols<\/span>))<\/span>\n    return<\/span> np<\/span>.<\/span>allclose<\/span>(<\/span>expected<\/span>,<\/span> actual<\/span>)<\/span>\n\ntestcases<\/span> =<\/span> [([<\/span>0<\/span>,<\/span> 2<\/span>-.<\/span>1<\/span>,<\/span> 2<\/span>+.<\/span>1<\/span>,<\/span> 5<\/span>],<\/span> [<\/span>1<\/span>,<\/span> .<\/span>1<\/span>,<\/span> .<\/span>1<\/span>,<\/span> 1<\/span>]),<\/span>\n             ([<\/span>0<\/span>,<\/span> 2<\/span>-<\/span>10e-8<\/span>,<\/span> 2<\/span>+<\/span>10e-8<\/span>,<\/span> 5<\/span>],<\/span> [<\/span>1<\/span>,<\/span> 10e-8<\/span>,<\/span> 10e-8<\/span>,<\/span> 1<\/span>]),<\/span>\n             ([<\/span>0.0<\/span>,<\/span>5<\/span>],<\/span> [<\/span>7.0<\/span>,<\/span>2<\/span>]),<\/span>\n             ([<\/span>0.0<\/span>,<\/span> 10<\/span>],<\/span> [<\/span>1e-4<\/span>,<\/span> 1e-4<\/span>]),<\/span>\n             ([<\/span>0.0<\/span>,<\/span> 1e5<\/span>],[<\/span>1.0<\/span>,<\/span> 4e-9<\/span>])]<\/span>\n\ndef<\/span> run_straight_tests<\/span>():<\/span>\n    print<\/span> "all pass?"<\/span>,<\/span> all<\/span>(<\/span>run_for<\/span>(<\/span>np<\/span>.<\/span>array<\/span>(<\/span>d<\/span>),<\/span> np<\/span>.<\/span>array<\/span>(<\/span>z<\/span>))<\/span>\n\nrun_straight_tests<\/span>()<\/span>\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n
        \n
        <\/div>\n
        \n
        all pass? True\n<\/pre>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        References<\/h3>\n<\/div>\n<\/div>\n<\/div>\n
        \n
        \n<\/div>\n
        \n
        \n
        The main references for this post come from:<\/p>\n
          \n
        1. Arbenz, Peter. Lecture Notes for Solving Large Scale Eigenvalue Problems, Chapter 5: Cuppen’s Divide and Conquer Algorithm.<\/a> The single most valuable piece of this is working Matlab code for solving the secular equation. I started with this code in octave, translated it to Python+NumPy, factored out some common pieces, renamed some variables, incorporated a different stopping condition, and massaged the variable names to match the discussion in this post.<\/li>\n
        2. Li, Ren-Cang. Solving Secular Equations Stably and Effciently<\/a>. This describes the theory behind LAPACK’s solution of the secular equation for the symmetric, tridiagonal eigenvalue problem. Arbenz follow the presentation of Li fairly closely. The LAPACK implementation – basically the LAPACK slaed<\/a> series of routines – is described by Rutter in a separate LAPACK Working Note<\/a>.<\/li>\n
        3. Demmel, James. Applied Numerical Linear Algebra, Section 5.3.3: Divide-and-Conquer. Great discussion of the connection between classical Newton’s method and the non-linear version we use. Great mathematical presentation of the setup for the approximation, but he punts on the details of solving the quadratic form.<\/li>\n<\/ol>\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":"
          It’s been way, way, way too long since I’ve posted. I haven’t been slacking though, I’ve merely been busy. Really. I decided to dive back into the SVD problem and look at an alternative to the QR based SVD computations. Namely, I’m going to give a breakdown of the divide-and-conquer approach to computing the SVD. […]<\/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-456","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\/456","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=456"}],"version-history":[{"count":1,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/456\/revisions"}],"predecessor-version":[{"id":457,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/posts\/456\/revisions\/457"}],"wp:attachment":[{"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/media?parent=456"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/categories?post=456"},{"taxonomy":"post_tag","embeddable":true,"href":"http:\/\/drsfenner.org\/blog\/wp-json\/wp\/v2\/tags?post=456"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}