[SciPy-user] do complex numbers default to double precision?

[Robert Kern]

Ryan Krauss wrote:
> Thanks to Nils and Robert for their quick responses.  This is definitely 
> one thing to love about SciPy.  I posted my null space question 
> yesterday and Robert responded in 20 minutes and Nils in 30.  I posted 
> my optimization question this morning and Nils reponded in 8 minutes!  
> It is almost like chatting with tech support. 

More pleasant than that, I hope! We have better "On Hold" music, too.

> You make SciPy a pleasure 
> to use!  Thanks.
> 
> Coincidentally, these two problems are linked and I don't know if my 
> numerical error problems are from fmin or the null space/svd stuff.  
> Once I find the input that drives my matrix to have a null space, I find 
> the vector that corresponds to the null space (assuming the sub-matrix 
> is only rank deficient by 1).  Then I combined the null space vector 
> with zeros that correspond with the boundary condition on one end of the 
> problem.  So, a 4x1 null space vector would give me an 8x1 full vector.  
> I then take the 8x8 full matrix which has a 4x4 submatrix whose 
> determinant is roughly 0 and multiply it by the 8x1 vector of the 
> boudnary conditions I just solved for.  This should then give me an 8x1 
> vector of the boundary conditions on the other end.  The problem is that 
> there are some elements of this second vector that should be 0 because 
> of the boundary conditions and they are actually of order 1e-10, if the 
> vector is normalized so that its magnitude is 1.  This physically means 
> that I have a cantilever beam with a free end that has just a little bit 
> of force and moment at the free tip.

A better function to minimize might be linalg.svdvals(A)[-1] (singular 
values are always sorted largest-to-smallest with LAPACK).

If the size of the elements of A are about 1e6 or so, the smallest 
singular value may be about 1e-10 and the result of multiplying it with 
a calculated null vector will have elements the size of 1e-10.

Come to think about it, a better nullspace function might be as follows:

def null(A, eps=1e-16):
     u, s, vh = scipy.linalg.svd(A)
     mask = (s/s[0] <= eps)
     rowspace = scipy.compress(mask, vh, axis=0)
     return scipy.conjugate(scipy.transpose(rowspace))

and a corresponding function to minimize:

def f(params):
     A = ...
     s = scipy.linalg.svdvals(A)
     return s[-1]/s[0]

That should be <~ 1e-16 regardless of the size of the matrix. Be sure to 
check the output of the minimizer; it may be stuck in a local minimum. 
Do some plots of f(params) to see how nasty your functions are.

You'll still get the small but > 1e-16 numbers when you multiply 
through, but if you've satisfied yourself that they ought to be zero, 
you can artificially set them to zero if you need to propagate the 
results through to other calculations.

Ain't floating point great?

-- 
Robert Kern
rkern at ucsd.edu

"In the fields of hell where the grass grows high
  Are the graves of dreams allowed to die."
   -- Richard Harter