Re: [boost] [Review][PQS] Review deadline

From: Leland Brown (lelandbrown_at_[hidden])
Date: 2006-06-09 02:50:32

John Phillips <phillips <at> delos.mps.ohio-state.edu> writes:

> The basis for chosing the units is
> instead attempting to maximize numeric stability, and so choosing units
> that keep things as close to order 1 as possible.

BTW, this post will be only of interest to those with a mathematics or 
numerical analysis background.  Others can safely skip it without missing 
anything.

This is a very interesting issue to me!  I'm not sure what you have in mind as 
far as keeping things close to order 1 improving numerical stability, but let 
me explain what comes to mind for me.  First, there may be a possibility of 
overflow or underflow with values significantly far from order 1, especially 
with single precision floating point values.  In my application, I use double 
precision for everything, so this doesn't concern me.

The second issue relates to things like inversion of ill-conditioned matrices, 
where a matrix may be nearly singular due to a large difference in the orders 
of magnitude of its eigenvalues.  Clearly this can cause severe numerical 
problems, and keeping things near order 1 helps prevent this.

In some cases FWIW I think this problem can be resolved in another way, or even 
seen as a symptom of a different problem, which can be addressed - as I'll 
explain.  And the discipline imposed by strongly-typed dimensional analysis may 
force the implementer of the algorithm to correct the problem.  (Of course, you 
can read this as "annoy the implementer because his mathematically correct 
algorithm won't compile," depending on your perspective :-) )

Consider, for example, a function that inverts a matrix A using an SVD, in 
order to avoid numerical problems.  Let's assume A is a symmetric, positive 
definite matrix composed of values having different physical dimensions - such 
as a covariance matrix of a set of parameters having different dimensions.  I 
don't use SVD (so forgive me if I display my ignorance); but as I understand, 
it would essentially find the eigenvalues and eigenvectors, and invert the 
eigenvalues; but for any very small eigenvalues it would replace the very large 
reciprocal value with a zero instead.  If we're solving Av=b, then loosly 
speaking the numerical problem is that many values of v solve the equation, 
some with very large coordinate values.  Geometrically, the SVD can be 
described as choosing from among these values of v the one closest to the 
origin.

What has always bugged me about this description is that "closest" requires a 
distance metric, like d=sqrt(x*x+y*y), but this is meaningless in a physical 
sense if the components x and y have different physical dimensions.  Sure, if 
you have numerical values in some particular units system, you can do the 
computation and get a number.  But what doesn't sit right with me is that if 
you change your units, you can get a different "closest" point!  So even 
without roundoff errors and the like, the answer that comes out of the 
algorithm will change depending on what physical units are chosen.  To me, it 
seems that a physical problem should have a particular solution regardless of 
the units chosen (aside from issues of computational accuracy).

Likewise, the eigenvalues of this matrix will have some sort of hybrid 
dimensions just like the distance metric above.  And the eigenvectors can vary 
based on choice of units, because even the notion of two vectors being 
orthogonal is dependent on the relationship among the units of the parameters 
defining the matrix.  Because of this, as you pointed out in another post, a 
strongly-typed dimensions library would not even be able to compile the SVD 
(unless all the elements of the matrix are of the same dimensionality).

And how do I set the threshold of what's a "very small" eigenvalue for the 
SVD?  Perhaps .000001?  In what units, since the matrix elements are of 
different dimensions?  Again, the meaning of my threshold changes depending on 
the relationship among my parameter units.

I used to think all this was just the price you pay for using an SVD.  But with 
the way I've constructed my matrix, it can be decomposed as:

   A = DCD

where D is a diagonal matrix of dimensioned quantities, and C is a 
dimensionless symmetric matrix.  It's easy, in fact, to find such a diagonal 
matrix D, with no a priori knowledge of the magnitude of the values, by taking 
the square roots of the diagonal elements of A.  Then C is just the correlation 
matrix corresponding to the covariance A, and its elements tend to be close to 
order 1.  Now the matrix to be inverted (C) is likely to avoid numerical 
problems caused by large differences in magnitude, and the solution will be 
independent of the units chosen - since they will just become multipliers on 
the elements of D.  The SVD/eigenvalue computation can be done without 
violating the strong typing, and the choice of units no longer impacts the 
numerical stability of the matrix inversion.

There may be a small penalty (such as calculating the square roots), but as 
compared to an SVD algorithm I think it's almost negligible.  I expect a 
similar technique could be found for cases of nondefinite matrices or non-
square matrices, and other types of matrix calculations.  In parts of the 
algorithm that are not "mixing" dimensions (as the distance formula above 
does) - where strongly-typed computations would compile - you're always adding 
apples to apples, so the choice of units won't influence numerical problems 
like summing values of very different magnitudes.  Most matrix multiplications 
should probably be type-safe, for instance, and thus independent of units as 
far as numerical stability.

These ideas would never have occurred to me when I was just manipulating 
matrices as set of numbers.  It wasn't until I started to think of them as 
strongly-typed physical quantities that this perspective emerged.  I think 
that's also part of the value of using a dimensional analysis library.

I hope what I've suggested here makes some sense.  There may be reasons it 
doesn't apply in certain situations.  But if anyone finds it helpful, it's 
something you can try.

Regards,

-- Leland