Re: Multivariate pdf of a normal distribution

[Mike Miller]

On Sun, 6 Nov 2005, Andy Adler wrote:

> On Sun, 6 Nov 2005, Mike Miller wrote:
>
> > When inv(A) exists,
> >
> > inv(A') = inv(A)'
> >
> > because the inverse equals the transpose of the adjoint divided by the
> > determinant.  The adjoint of the transpose equals the transpose of the
> > adjoint and the determinant is the same for A and A'.  The identity really
> > follows from the fact that the determinant is not affected by the
> > transpose operation -- the adjoint is a collection of determinants.
> >
> > I agree that this must be in most matrix algebra books.
>
> Yes, that is true in theory. But the error is to assume that all
> mathematical identities remain intact when done on a finite
> precision computer.
>
> Try it:
> A= rand(2); inv(A')- inv(A)'
> ans =
>
>    7.1054e-15   -7.1054e-15
>    8.8818e-16    0.0000e+00


That wasn't an error that anyone here was making.  The original question 
showed an example where there was a small difference.

Related to this, we also were discussing recently the problem with eigen 
decomposition of singular matrices.  For example:


octave:6> [E,L] = eig(ones(4))
E =

   -3.2026e-01    7.0711e-01    3.8397e-01    5.0000e-01
   -3.2026e-01   -7.0711e-01    3.8397e-01    5.0000e-01
   -2.2276e-01   -5.9234e-17   -8.3689e-01    5.0000e-01
    8.6328e-01   -4.4132e-19    6.8941e-02    5.0000e-01

L =

  -0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   4.00000

octave:7> A = E*sqrt(L)
A =

   0.00000 - 0.00000i   0.00000 + 0.00000i   0.00000 + 0.00000i   1.00000 + 
0.00000i
   0.00000 - 0.00000i  -0.00000 + 0.00000i   0.00000 + 0.00000i   1.00000 + 
0.00000i
   0.00000 - 0.00000i  -0.00000 + 0.00000i  -0.00000 + 0.00000i   1.00000 + 
0.00000i
   0.00000 + 0.00000i  -0.00000 + 0.00000i   0.00000 + 0.00000i   1.00000 + 
0.00000i

octave:8> A*A'
ans =

  1.0000  1.0000  1.0000  1.0000
  1.0000  1.0000  1.0000  1.0000
  1.0000  1.0000  1.0000  1.0000
  1.0000  1.0000  1.0000  1.0000

The problem here is that numerical imprecision causes A to be complex when 
we want it to be real.  If the matrix really is singular (as this in the 
example above), we can use real() to get a good answer:

octave:9> B=real(A)
B =

   0.00000   0.00000   0.00000   1.00000
   0.00000  -0.00000   0.00000   1.00000
   0.00000  -0.00000  -0.00000   1.00000
   0.00000  -0.00000   0.00000   1.00000

octave:10> B*B'
ans =

  1  1  1  1
  1  1  1  1
  1  1  1  1
  1  1  1  1

But for a function that always works, we'll have to set some tolerance 
limit and have the function fail when the smallest eigenvalue [min(L) 
above] is less than some value.  What is a good choice of a value?  Maybe 
-1.0e-13?  Anyone?  How far off can we be when the precise min(L) is zero?

Mike



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------