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

From: Leland Brown (lelandbrown_at_[hidden])
Date: 2006-06-09 01:28:24

Andy Little wrote (in an earlier post):

> OK, what I was going onto say was that this is only the beginning of the 
> complexities involved in strong type checking on physical quantities .

Yeah, and wait till you hear this one!  It gets worse - though I also have a 
solution that works for me.  Read on, if you dare :-)

Warning - this post is primarily for the more mathematically inclined - though 
I think it brings up an important issue that anyone wanting to implement a 
matrix library on dimensioned quantities needs to think about.

John Phillips <phillips <at> mps.ohio-state.edu> writes:
> 
>    Maybe I'm not getting it yet, but I don't think this gets you out of 
> the woods.
>    The problem I see is that matrices are used for many things other 
> than just transformations. For example, a user might want to diagonalize 
> a matrix, or find the eigenvalues, or any number of other things. If so, 
> having units on the matrix becomes problematic, I think.

With certain restrictions, this can be done.  But it brings up other problems.  
Consider:

It's a basic tenet of matrix algebra that for a square matrix A (as long as it 
has an inverse), both of these are true:

   A * inv(A) = I
   inv(A) * A = I

where I is the identity matrix of the same size as A.  Now suppose A is a 
matrix of elements having different dimensions.  For instance, if A is a 2x2 
covariance matrix of position and velocity values, it might have units like

   A:
      m^2     m^2/s
      m^2/s   m^2/s^2

Its inverse has well-defined dimensions, with the units you might expect, like:

   inv(A):
      1/m^2   s/m^2
      s/m^2   s^2/m^2

But if you multiply A * inv(A) or inv(A) * A, you'll get two different 
products, neither of which is the identity matrix!!!  The results are:

   A * inv(A) =
      1    0 s      (I call this "I1" below)
      0/s  1

   inv(A) * A =
      1    0/s      (I call this "I2" below)
      0 s  1

   I =
      1    0
      0    1

Of course, we all learned in high school physics that "0" and "0 seconds" are 
*not* the same thing, and Andy's library faithfully embodies that truth.  Thus, 
we have three *different* matrix values above.  To make matters worse - for the 
same reason, neither A * I = A nor I * A = A.  In fact, neither expression
(A * I) nor (I * A) will compile, because they both involve addition of 
quantities with different dimensions.  This addition would be (correctly, I 
think) not allowed by the library, even though in this particular case, one of 
the addends has a numerical value of zero.  To multiply A by the identity, we 
would need (I1 * A) or (A * I2).  What a mess!

BTW, I believe John's example of Guassian elimination may fail not only because 
of things like pivoting (as he points out), but even without that, because of 
this issue - addition or subtraction of values where one is numerically exactly 
zero but has the wrong dimensions.

So what's a programmer to do?  In my case, I was already using something akin 
to "t3_quantity" for these "mixed" matrices, and in that context I came up with 
the following solution:

I introduced a special value that a t3_quantity can hold, called an 
adimensional zero.  Normally, a value of my t3_quantity has specific dimensions 
(or is specifically dimensionless), even if its numerical value happens to be 
zero.  An adimensional zero is a special case.  What this means is that its 
dimensions are undefined, ambiguous.  Semantically, this means it can be added 
to a value of any dimension (zero + x = x, of course, with the dimensions of 
the result the same), and when multiplied by any value the result is another 
adimensional zero (zero * x = zero).

I define a constant called "zero," and then define the identity matrix as

   I =
       1    zero
      zero   1

Now (I1 == I) and (I2 == I) both evaluate to true, (A * I) and (I * A) both 
compile and give A as the result, etc.  Problem solved!  It works out quite 
well in all my linear algebra.  But sadly, at the price of increasing the run-
time overhead already required by my "t3_quantity."

Anybody got a better idea?

-- Leland