On Mon, Sep 7, 2015 at 5:37 PM, Carsten Knoll <[email protected]> wrote: > Hi, > > I want to write a routine to determine the rank of a symbolic matrix and > it seems to be a quite tricky task. > > The routine is required to be applicable for matrices whose entries are > "big expressions", hence calling simplify is not an option. > > > My strategy: extracting all occurring symbols and replacing them by > random numbers. This way, I arrive at a numeric matrix, which might > contain very small numbers which are nevertheless different from zero. > Then, I call evalf with different precisions and look at the singular > values. My assumption was, that the singular values which are > "practically 0" must become smaller with rising precision while singular > values whose true value is small but greater zero wont change much.
I don't know if that is going to be true. I believe what you are dealing with is related the so-called table-maker's dilemma https://en.wikipedia.org/wiki/Rounding#Table-maker.27s_dilemma. There's no way to tell how many digits you need to get a "zero" value below a certain epsilon. That isn't to say that there aren't good heuristics for doing this, some of which are implemented in SymPy (see equal() for instance). Aaron Meurer > > I tried three ways to calculate the singular values, all of which lead > to problems: > > 1. mpmath: I could not figure out how to increase precision > > > 2. numpy: precision is limited by data type (np.float64) > > > 3. direct calculation by (M.T*M).eigenvals() > This often works as expected but fails for the following expample: > > > M = Matrix([ > [1, 0, 0], > [1, sin(3)**50, 1], > [0, 0, 0]]) > > due to the 0-row, one singular value should be exactly zero. > > But I get the following: > > In [14]: Matrix((M.T*M).evalf(prec=20).eigenvals().keys()).evalf(20) > Out[14]: > Matrix([ > [ 2.6180339887498948482 - 0.e-27*I], > [ 0.3819660112501051518 - 0.e-28*I], > [2.4515092575814173921e-102 + 0.e-125*I]]) > > In [15]: Matrix((M.T*M).evalf(prec=50).eigenvals().keys()).evalf(50) > Out[15]: > Matrix([ > [ 2.6180339887498948482045868343656381177203091798058 - 0.e-56*I], > [ 0.38196601125010515179541316563436188227969082019424 - 0.e-58*I], > [2.4515092575814173921301544696749858514095930873163e-102 + 0.e-156*I]]) > > > This contradicts my assumption that with increasing precision the last > singular value should get smaller (i.e. closer to 0). > > Interestingly, if I swap M and M.T (which leaves the singular values of > a quadratic matrix invariant by definition ) I get the expected result: > > > In [16]: sp.Matrix((M*M.T).evalf(prec=20).eigenvals().keys()).evalf(20) > Out[16]: > Matrix([ > [ 0], > [2.6180339887498948482], > [0.3819660112501051518]]) > > > Questions: > > 1. Is this a viable strategy for rank calculation? Are there preferable > alternatives? > > 2. Is it possible to use mpmath to calculate arbitrary precise singular > values? > > 3. If not, how could I alternatively calculate the eigenvalues of > (M.T*M) with arbitrary precision? > > > 4. Is the non-convergence to zero for the last singular value in the > above example a bug or is my expection wrong? > > > Thanks in advance and best regards, > Carsten. > > > -- > You received this message because you are subscribed to the Google Groups > "sympy" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to [email protected]. > To post to this group, send email to [email protected]. > Visit this group at http://groups.google.com/group/sympy. > To view this discussion on the web visit > https://groups.google.com/d/msgid/sympy/55EE11C6.8090309%40gmx.de. > For more options, visit https://groups.google.com/d/optout. -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at http://groups.google.com/group/sympy. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6%2BZayBZ_9m%2B%2BO5FN3GyhuRzNgzF5AtNnYjCcHnsp%2B%3D3qA%40mail.gmail.com. For more options, visit https://groups.google.com/d/optout.
