7
$\begingroup$

I have an application in which I am computing a quantity which is approximated by an average over $M$ points. In theory, the average converges to the correct quantity when $M$ is infinite. In practice, the computation involves a $M\times M$ symmetric positive definite matrix, whose condition number gets huge when $M$ increases. I compute the condition number as the ratio of the biggest and smallest eigenvalues given by a selfadjoint eigensolver.

Therefore, the average I'm computing does not seem to reach stationarity when I increase $M$, but instead starts deviating quite a bit from what I feel is the correct value when $M$ is really big.

How can I choose the best $M$? Is there a way to put some error bars on the result I get, based on the condition number? What is a bad condition number for the $LDL^T$ cholesky decomposition and for a SPD matrix, e.g. when should I start worrying?

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ So...why are you using a self-adjoint eigensolver and not a singular value decomposition? Also, what is the linear algebra problem you're trying to set up? $\endgroup$
    – Geoff Oxberry
    Jan 21, 2012 at 20:30
  • $\begingroup$ @Geoff Oxberry: His matrix is SPD, so the eigenvalues equal the singular values. In fact, the self adjoint eigensolve should be significantly faster. $\endgroup$
    – Jack Poulson
    Jan 21, 2012 at 21:37
  • $\begingroup$ @JackPoulson: That was where I was going with the first point, but I'm surprised the self-adjoint eigensolve is faster. I suppose for the condition number you only need two eigenvalues, and extremal eigenvalues are easier to calculate. $\endgroup$
    – Geoff Oxberry
    Jan 21, 2012 at 21:47
  • 2
    $\begingroup$ @GeoffOxberry: It is essentially a result of SVD implementations almost always not exploiting symmetry. I would implement self adjoint SVD as a wrapper around a self adjoint eigensolver (even if not positive definite). $\endgroup$
    – Jack Poulson
    Jan 21, 2012 at 22:07
  • $\begingroup$ @GeoffOxberry I am writing a multivariate normal distribution, so the condition number is just for debug purposes. All I need is determinant and solving for Ax=b, therefore I feel cholesky is perfectly suited to my needs. $\endgroup$
    – yannick
    Jan 22, 2012 at 16:42

4 Answers 4

Reset to default
7
$\begingroup$

I'm confused as to why you're saying the $LDL^T$ decomposition gives you your eigenvalues; one can certainly compute the determinant (and inertia) of a symmetric matrix cheaply from it, but it does not directly yield the eigenvalues themselves, as $L \cdot L^T$ is not a similarity transformation unless $L=I$. Moreover, I would suggest using the Cholesky (i.e., $LL^T$) decomposition for SPD matrices instead of $LDL^T$.

If it does turn out that the condition number really is terrible, then you are still in decent shape, as preconditioning SPD matrices is usually not very challenging, especially for diagonally-dominant matrices.

EDIT: Since you updated your question, I will update my response.

If you condition number is really getting as high as, say, $10^{10}$ (so that you will lose roughly 10 decimal digits, as @aeismail mentioned), for matrices which you can compute the dense eigensolution on a regular machine, then I would start by making sure that no numerical atrocities are being committed. In particular, I would make sure that your self-adjoint matrix is not being formed as the symmetric product of a matrix with itself, e.g., $A := B B'$, as $A$'s condition number will be the square of $B$'s condition number. It would be much better to work directly with $B$.

On the other hand, there are often ways to cheaply improve the quality of your solution. For instance, if you were solving a linear system with your $LDL^T$ decomposition, you could "iteratively refine" your solution by running a handful of iterations of GMRES (or your favorite iterative method) with your solution from the direct solve as the initial guess. You haven't stated what you want to do with your system, so I can't say much more.

Improve this answer
$\endgroup$
4
  • $\begingroup$ On the other hand, one might want to stick with the $\mathbf L\mathbf D\mathbf L^\top$ decomposition if the square roots are found to be expensive... $\endgroup$
    – J. M.
    Jan 21, 2012 at 6:27
  • $\begingroup$ @J.M. Square root costs about the same as division on most hardware (akuvian.org/src/mubench_results.txt, agner.org/optimize/instruction_tables.pdf). And for dense(r) matrices, the diagonal operation is a lower order term that becomes insignificant even at very small sizes. $\endgroup$
    – Jed Brown
    Jan 21, 2012 at 14:30
  • 1
    $\begingroup$ @J.M. The situation is slightly more subtle: implementations of $LDL^T$ typically pivot, whereas it is rare for a Cholesky decomposition to pivot. Also, even though an unpivoted $LDL^T$ factorization almost performs the same operations as an unpivoted Cholesky factorization, the outer-product updates in a right-looking Cholesky map nicely to the BLAS routine dsyrk/zherk, whereas there is no such routine for the symmetric update in right-looking $LDL^T$. $\endgroup$
    – Jack Poulson
    Jan 21, 2012 at 17:52
  • $\begingroup$ I agree with @Jack Poulson, see [link]eigen.tuxfamily.org/api/classEigen_1_1LDLT.html[link] which is the library I use. However I assumed that the diagonal elements of D are the eigenvalues. Thanks for the hint. Correction above. $\endgroup$
    – yannick
    Jan 21, 2012 at 19:27
4
$\begingroup$

To comment on the question about the error bounds for the matrix: it's a little tricky to specify strict bounds, but there is a heuristic rule in numerical linear algebra that a condition number of $M$ will lead to a loss of approximately $\log_b M$ base-$b$ digits of precision in the results of your calculation involving the matrix $M$. So, for instance, a condition number $M = 100$ means you would lose $\log_{10}(100) = 2$ digits of precision.

Improve this answer
$\endgroup$
2
  • $\begingroup$ You need to be more precise: for condition number $\kappa$, one (heuristically) loses $\lfloor\log_b \kappa \rfloor$ or so base-$b$ digits. $\endgroup$
    – J. M.
    Jan 22, 2012 at 16:39
  • $\begingroup$ You have provided a link to a textbook. On which page of the textbook is the thing you're trying to show us? Also, due this heuristic rule apply to any matrix computation involving $M$? What about multiplication by the identity? $\endgroup$
    – Nike Dattani
    Apr 13, 2021 at 21:23
2
$\begingroup$

There is also a brute-force way to cope with ill-conditioning for not very large matrices: you can try to increase the precision of the machine arithmetics. If you use C or other supported language the obvious choice is GNU GMP library. But it can be made much simpler (and more computationally expensive) in Mathematica using SetPrecision[...] function.

Of course, the more the condition number the more digits of precision you'll need to get the satisfactory result.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ A way to use this strategy without going all-out is to use some form of iterative refinement, where your residuals are calculated in the extended arithmetic and everything else done in double. I have performed calculations with horribly conditioned matrices using this technique before to acceptable results. Since extended precisions are not supported natively on almost any computer hardware, you usually come out winning this way (a few solves in double will be faster than one solve in extended precision). $\endgroup$
    – Reid.Atcheson
    Jan 24, 2012 at 14:07
2
$\begingroup$

For rigorous, but probably weak error bounds, you can consult Nick Higham's Accuracy and Stability of Numerical Algorithms, 2nd edition. A posteriori error estimates are probably only good to first order, but might be useful. Don Estep does work on these mostly with regards to solutions of partial differential equations, but he does give a linear algebra example buried in slide 72 of this presentation, and a much more fleshed out example is given in slides 53-56 of this presentation. These estimates can do a pretty good job if the error is small, but tend not to be so accurate for larger errors.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.