How many digits of of accuracy do I expect the have the solution $x$ of $||Ax-b||=0$

condition numberfloating pointleast squaresmatrix decompositionnumerical linear algebra

A least-square problem $||Ax-b||=0$ is solved using a backward stable algorithm (In my case, QR decomposition using householder projectors). The condition number is $\kappa(A)=10^5$.

If the problem is solved using double-precision floating-point arithmetic($10^{-16}$), how many digits of accuracy I should expect the solution $x$ to have?

What I have tried:

I have tried to use this inequalities form the trefethen book page 131-Numerical Linear Algebra

And also I used this table from the same book.

I also found that the exponent of $\kappa(A)=10^5$ means that I will lose 5 digits of accuracy. I would conclude that the accuracy is $10^{-11}$

My doubt: is it this reasoning is ok?. And the second question is:
When do I have to take into account $\kappa(A)^{2}$ rather than only $\kappa(A)$?

Best Answer

To answer your first question, a rule of thumb for a linear system of equations is that if the unit roundoff is $10^{-a}$ and the condition number is $10^b$, then you can expect about $a-b$ digits of accuracy in your answer. This is because, for a backwards stable algorithm, the numerical errors made during your computation can be viewed as a relative perturbation of size $\approx 10^{-a}$ of the initial data. Then the forward error is smaller than the backwards error times condition number so the forward relative error for $x$ is $\approx 10^{b-a}$, leading to $\approx a-b$ correct digits.

A least squares problem has additional subtleties over a pure linear system. For a formal derivation, Trefethen and Bau explain how the $(\kappa(A))^2$ appears in Eqs. (18.13-18.16) of the book you quote. Is their explanation good enough or do you have additional questions?

Regarding the $(\kappa(A))^2$ term in the error bound, note that the $(\kappa(A))^2$ term is tamped by two additional factors. if $\eta$ is on the order of $\kappa(A)$ (which will "usually" happen if $y$ is "randomly chosen"), then $(\kappa(A))^2 / \eta \approx \kappa(A)$. Additionally, if the angle between $y$ and $b$ is small (that is, if $Ax$ is a good approximation for $b$), then $\tan\theta$ will be $\approx 0$ and $(\kappa(A))^2 \tan(\theta)/ \eta$ will be small or at least on the order of $\kappa(A)$. Thus, the $(\kappa(A))^2$ term in the error bound is usually not as bad as it seems, and the rule of the thumb from the first paragraph is usually valid for least squares problems as well.

For some vague intuition as to why we might not be surprised $(\kappa(A))^2$ pops up, remember that the least squares problem is mathematically equivalent to the normal equations $A^\top A x = A^\top b$, which has condition number $\kappa(A^\top A) = (\kappa(A))^2$. Normally, the least-squares problem is significantly better conditioned than the normal equations, but in the worst case (which as established in the previous paragraph is a somewhat special situation) the conditioning becomes the same.

Related Solutions

[Math] Why do we say SVD can handle singular matrix in least-squares? Comparison of SVD and QR decompositions

SVD can handle rank-deficiency. In your example, there was a bug. The diagonal matrix D has a near-zero element and you need use pseudoinverse for SVD, i.e. set the 2nd element of 1./diag(D) to 0 other than the original huge value (10^14). You should find SVD and QR have equally good accuracy in this case. For more information, see this document http://www.cs.princeton.edu/courses/archive/fall11/cos323/notes/cos323_f11_lecture09_svd.pdf

Is There Any Good Way of Checking Floats for Approximate Equality

The easiest way to check for float near-equality which scales is something like $$ {|a-b|\over|a|+|b|}<\epsilon $$ where $\epsilon$ (epsilon) is the tolerance.

In code that would be

package main

import (
    "fmt"
    "math"
)

func within(a, b, tol float64) bool {
    d := math.Abs(a - b)
    if d == 0 {
        return true
    }
    // Since d =/= 0 we know a, b =/= 0 so division by zero will not happen
    return d/(math.Abs(a)+math.Abs(b)) < tol
}

func main() {
    checks := [][3]float64{
        {0, 0, 0.1},
        {100000, 100001, 0.0001},
        {0.00001, 0.000014, 0.0001},
        {0.00001, 0.00001000001, 0.0001},
    }
    for _, c := range checks {
        fmt.Printf("%g %g %g %t\n", c[0], c[1], c[2], within(c[0], c[1], c[2]))
    }
}

The output looks like this:

0 0 0.1 true
100000 100001 0.0001 true
1e-05 1.4e-05 0.0001 false
1e-05 1.000001e-05 0.0001 true

Run online here.

A caveat with this method is that it doesn't work for the case that you expect that the output will be zero, for example if you expect $a=0$ then the above calculation doesn't give useful results. In that case you need to have an absolute value for the tolerance, $$|b|<\epsilon,$$ instead of the above.

Best Answer

Related Solutions

[Math] Why do we say SVD can handle singular matrix in least-squares? Comparison of SVD and QR decompositions

Is There Any Good Way of Checking Floats for Approximate Equality

Related Question