[SciPy-user] do complex numbers default to double precision?
Ryan Krauss
ryanfedora at comcast.net
Thu Jun 30 11:36:36 CDT 2005
So what should a good plot of the function I want to minimize look
like? To my untrainded eye, the determinant seems like it is slightly
better conditioned as magnitude of the eigenvalue increases, while the
ratio of the minimum to the maximum svd looks better for smaller
magnitudes. The svd approach is consdierably faster for the same
convergence criteria - mainly because the determinant based method
spends a long time looking for the second eigenvalue.
Ryan
Robert Kern wrote:
> 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?
>
More information about the SciPy-user
mailing list